1 Rプログラミングの基礎

1.1 RStudioの基本

以下のサイトからRとRStudioのインストーラをダウンロードし、インストール
R: https://cran.ism.ac.jp/ または https://cran.r-project.org/mirrors.html
RStudio: https://posit.co/download/rstudio-desktop/

左上のSource paneでR scriptにコードを記述する。
コードを選択してCmd+Return またはOption+Return でコードが実行される。
Cmd+Option+P で直前に実行したコードを再実行できる、コードを一部変更して確認する際に便利。

右上のEnvironment paneには作成したオブジェクトが表示される。オブジェクトの中身を表示させたり、オブジェクトの削除も可能。
右下のFile paneではフォルダ内のファイル、作成したグラフ、機能の説明などが表示される。Packagesパネルではパッケージのインストールとロードも可能。
基本的に1件のデータ解析ごとに専用のフォルダを作成し、それをR projectファイルと関連させて管理する。

1.2 使用するパッケージのロード

library(pacman)
p_load(tidyverse, readxl, openxlsx)  # 一度に複数のパッケージをロードできる

1.3 Rにおける計算の基本

1.3.1 基本演算

Rは動的プログラミング言語（インタプリタ型言語）であり、Rコードを実行するためにコンパイルする必要は無く、1行ずつコードが解釈されて実行される。
四則演算は+, -, *, /で可能。ハッシュタグ #はコードにおけるコメントを表す。

1 + 2  # #記号より右は無視される
## [1] 3
4 - 1
## [1] 3
12 * 6
## [1] 72
5 ^ 3 # べき乗
## [1] 125
20 / 3
## [1] 6.666667
20 %/% 3 # 整数商
## [1] 6
20 %% 3  # 剰余（mod）
## [1] 2

コロン:は1刻みで連続な数を簡単に作成できる。他の四則演算と組み合わせることも可能。

1:5
## [1] 1 2 3 4 5
-2:3
## [1] -2 -1  0  1  2  3
-1:-10
##  [1]  -1  -2  -3  -4  -5  -6  -7  -8  -9 -10
1:5 + 2
## [1] 3 4 5 6 7
1:5 - 2
## [1] -1  0  1  2  3
1:5 / 2
## [1] 0.5 1.0 1.5 2.0 2.5
1:5 * 2
## [1]  2  4  6  8 10
1:5 * 1:5
## [1]  1  4  9 16 25

1.3.2 関数

Rには基本的な数学操作を行う関数や、データの統計的情報を知るための集計関数が最初から用意されている。
関数に渡すデータは引数（ひきすう）と呼ばれる。データそのもの、Rオブジェクト、他の関数の結果、いずれも引数になり得る。 下記のround関数におけるdigitsのように、関数で指定する対象の引数を仮引数（parameter）と呼び、仮引数に実際に 渡される値のことを実引数（argument）と呼ぶ。誤解が生じない場合は2つを区別せずに引数とする。

round(3.1415, digits = 2)  # 数値を丸める、digitsは小数点以下の桁数
## [1] 3.14
factorial(4)  # 階乗を計算（4*3*2*1）
## [1] 24
factorial(round(3.14, digits = 0)) # 関数の結果が引数になる場合、最も内側の関数から外側に向かって関数が実行される。
## [1] 6
x <- 1:50
length(1:20) # データの個数
## [1] 20
max(x) # 最大値
## [1] 50
which.max(x) # 最大値の位置（インデックス）
## [1] 50
min(x) # 最小値
## [1] 1
which.min(x) # 最小値の位置
## [1] 1
mean(x) # 平均値
## [1] 25.5
median(x) # 中央値
## [1] 25.5
sd(x) # 標準偏差
## [1] 14.57738
var(x) # 分散
## [1] 212.5
range(x) # 範囲
## [1]  1 50
quantile(x) # 分位点
##    0%   25%   50%   75%  100% 
##  1.00 13.25 25.50 37.75 50.00
summary(x) # 要約統計量
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   13.25   25.50   25.50   37.75   50.00

関数の引数の指定は省略可能。またargs()関数で関数の引数名とデフォルト値を調べることができる。

round(3.1415, 2)
## [1] 3.14
args(round)
## function (x, digits = 0) 
## NULL
# 数字の桁数を揃える
pi
## [1] 3.141593
format(pi, nsmall = 8) # 少数点以下8で揃える
## [1] "3.14159265"
format(pi, digits = 2) # 数字2桁以下を切り捨てる
## [1] "3.1"
# 0埋め（zero-padding）
formatC(5:12, width = 2, flag = "0")
## [1] "05" "06" "07" "08" "09" "10" "11" "12"

連続な数を生成する関数

seq(1, 10, length = 5)  #1〜10を5分割で
## [1]  1.00  3.25  5.50  7.75 10.00
seq(1, 10, by = 2)  #1から2ずつ増やして10まで
## [1] 1 3 5 7 9
seq(as.Date("2016-01-01"), as.Date("2016-01-05"), by = "day")  # 日付も生成できる
## [1] "2016-01-01" "2016-01-02" "2016-01-03" "2016-01-04" "2016-01-05"
rep(1:3, times = 2) # (1, 2, 3)を2回
## [1] 1 2 3 1 2 3
rep(1:3, length = 5) #(1, 2, 3)を5個になるまで繰り返す
## [1] 1 2 3 1 2
rep(1:3, each = 3) #(1, 2, 3)の要素を3個ずつ
## [1] 1 1 1 2 2 2 3 3 3
replicate(10, 1+1) ## replicateはコマンドを複数回実行して結果をベクトルとして保存
##  [1] 2 2 2 2 2 2 2 2 2 2

無作為標本抽出のための関数

sample(x = 1:6, size = 4)  # ベクトルxからsize個の要素を取り出す
## [1] 2 4 3 5
sample(x = 1:6, size = 4, replace = TRUE)  # replace = TRUEは取り出した要素を元に戻す（要素の重複が許される）
## [1] 3 6 4 6
x <- rnorm(100, mean = 0, sd = 1) # 平均0, 標準偏差1の正規分布からランダムに値を100個抽出

割合を指定して無作為抽出するにはsample_frac()を使う。ただし入力はデータフレーム。

sample_frac(data.frame(x = 1:20), size = 0.2)  # 1〜20から20%の要素を取り出す
##    x
## 1  9
## 2 19
## 3 15
## 4  7

1.3.3 オブジェクト

Rではオブジェクトの中にデータを格納して保存することができる。オブジェクトを作るには名前を選び、 代入演算子 <-を使用する。 使用中のオブジェクトはEnvironmental paneに表示される。同じ名前のオブジェクトは上書きされる。

a <- 1 # RStudioではOption + (-)キーが代入演算子のショートカット
3 -> b # 向きが逆でも大丈夫
a
## [1] 1
b
## [1] 3
ls() # すでに使用しているオブジェクト名がわかる
## [1] "a" "b" "x"
assign("C", 7, envir = globalenv()) # assign関数は代入演算子 -> と同じ、ただし環境を指定できる

代入操作をカッコで括ると代入とコンソールでの表示を同時に行う。

(c <- 1 + 3)
## [1] 4

1.3.4 データ型

Rのオブジェクトにはbasic type, mode, class,の3つの型がある：

  • basic typeはデータを保存できるオブジェクトの型であり、vector, list, function, language, expression, environment等がある。 basic typeを調べるには関数typeof()を使用する。typeof()はvector objectに対して使用するとmode型を返す。
  • modeはオブジェクトに格納された要素に対する型であり、logical, numeric, complex, character, rawがある。numericはさらにdoubleとintegerに分かれる。 modeを調べたり、変更するには関数mode()を使用する。
  • classはオブジェクトの属性（attributes）に対する型であり、matrix, array, factor, data.frame等がある。 classが定義されていない場合は、basic typeやmodeの型を継承する。 classを調べたり、変更するには関数class()を使用する。

1.3.5 ベクトル

Rでは基本的にデータをベクトル単位で扱う。データをc()で囲むとベクトルオブジェクトになる。ベクトルの要素単位で計算を実行する。

v1 <- c(8,7,2,5,1,3)
v2 <- c(1:6)
v1
## [1] 8 7 2 5 1 3
v2
## [1] 1 2 3 4 5 6
v1 - 2
## [1]  6  5  0  3 -1  1
v1 / 2
## [1] 4.0 3.5 1.0 2.5 0.5 1.5
v1 - v2
## [1]  7  5 -1  1 -4 -3
v1 * v2
## [1]  8 14  6 20  5 18

ベクトルは1つの同じ型のデータしか格納できない。

int <- c(1L, 5L)  # 数字の後にLを付けると整数integerになる
text <- c("ace", "hearts")
int
## [1] 1 5
typeof(int)
## [1] "integer"
text
## [1] "ace"    "hearts"
typeof(text)
## [1] "character"

1.3.6 行列

行列は数学の行列とほぼ同じ、行と列からなる2次元の数値の配列。

die <- 1:6
m <- matrix(die, nrow = 2) ## 行列を作成、次元属性が付与される
m
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
m <- matrix(die, nrow = 2, byrow = T)
m
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
dimnames(m) <- list(c("r1","r2"), c("c1","c2","c3")) ## dimnames属性で行列の行と列に名前をつけることもできる
m
##    c1 c2 c3
## r1  1  2  3
## r2  4  5  6

1.3.7 配列

配列（array）は行列の拡張版であり、行列は配列の特殊な形。
配列は同じサイズの行列を数枚重ねたものであり、行と列以外に「層」を持つ。

1.3.7.1 数列から直接、配列を作成

ar <- array(c(11:14, 21:24, 31:34), dim = c(2,2,3))  ## array関数はn次元配列を作成, dim属性を指定する
ar
## , , 1
## 
##      [,1] [,2]
## [1,]   11   13
## [2,]   12   14
## 
## , , 2
## 
##      [,1] [,2]
## [1,]   21   23
## [2,]   22   24
## 
## , , 3
## 
##      [,1] [,2]
## [1,]   31   33
## [2,]   32   34

1.3.7.2 行列を重ねて配列を作成

3行4列の行列を4層重ねて配列を作成。

Mat1 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat2 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat3 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat4 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Array1 <- array(c(Mat1, Mat2, Mat3, Mat4), dim = c(3, 4, 4))
class(Array1)
## [1] "array"
Array1
## , , 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]    2    5    7   11
## [2,]    2    5    6    6
## [3,]    4    8   11    9
## 
## , , 2
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    7    2   11
## [2,]    3   10    3    4
## [3,]    2    7    4    1
## 
## , , 3
## 
##      [,1] [,2] [,3] [,4]
## [1,]   10    5    5    4
## [2,]   10    6   10   11
## [3,]    6   11    3   11
## 
## , , 4
## 
##      [,1] [,2] [,3] [,4]
## [1,]    7    7   10   10
## [2,]    4    8    8   10
## [3,]    1    4    4    1

配列名[行番号, 列番号, 層番号]で要素を指定。

Array1[3, 1, 2]  # 2番目の行列の3行1列目の要素を抽出
## [1] 2
Array1[ , , 3]  # 3番目の行列を抽出
##      [,1] [,2] [,3] [,4]
## [1,]   10    5    5    4
## [2,]   10    6   10   11
## [3,]    6   11    3   11
Array1[1:2, 1:2, ]  # 各層から1-2行目と1-2列目の要素からなる2×2行列を抽出
## , , 1
## 
##      [,1] [,2]
## [1,]    2    5
## [2,]    2    5
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    1    7
## [2,]    3   10
## 
## , , 3
## 
##      [,1] [,2]
## [1,]   10    5
## [2,]   10    6
## 
## , , 4
## 
##      [,1] [,2]
## [1,]    7    7
## [2,]    4    8
Array1[2, , ]  # 全ての行列の2行目を抽出、結果は行列になることに注意
##      [,1] [,2] [,3] [,4]
## [1,]    2    3   10    4
## [2,]    5   10    6    8
## [3,]    6    3   10    8
## [4,]    6    4   11   10

1.3.8 リスト

リストは任意のRオブジェクトを1次元のグループにまとめる。

list1 <- list(100:130, "R", list(T, F)) 
names(list1) <- c("sequence", "character", "logical")  ##  namesで名前属性をつけることができる
list1
## $sequence
##  [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
## 
## $character
## [1] "R"
## 
## $logical
## $logical[[1]]
## [1] TRUE
## 
## $logical[[2]]
## [1] FALSE
str(list1)  # まとめているオブジェクトの型はstr()関数で見ることができる
## List of 3
##  $ sequence : int [1:31] 100 101 102 103 104 105 106 107 108 109 ...
##  $ character: chr "R"
##  $ logical  :List of 2
##   ..$ : logi TRUE
##   ..$ : logi FALSE
list1[1] # 番号で要素を選択、リスト形式で返される
## $sequence
##  [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
list_len3 <- vector("list", length = 3)  # vector()で任意の長さのリストを作成できる

1.3.9 データフレーム

データフレームはリストの2次元バージョンであり、R版のExcelスプレッドシート。データフレームはベクトルをまとめて2次元の表にし、個々のベクトルは列になる。 それぞれのベクトルは異なる型のデータを格納できるが、同じ長さにする必要がある。言い換えるとデータフレームとは長さの等しいベクトルの名前付きリスト。データ分析のためには最も便利なストレージ構造。

df <- data.frame(face=c("ace","two","six"), suit=c("clubs","hearts","diamonds"),  # names, class, row.names属性をもつ
                 value=c(1,2,3), stringsAsFactors = F)
df
##   face     suit value
## 1  ace    clubs     1
## 2  two   hearts     2
## 3  six diamonds     3
str(df)
## 'data.frame':    3 obs. of  3 variables:
##  $ face : chr  "ace" "two" "six"
##  $ suit : chr  "clubs" "hearts" "diamonds"
##  $ value: num  1 2 3
typeof(df)
## [1] "list"
class(df)
## [1] "data.frame"
# dput()関数はオブジェクトをテキスト形式で出力する。defaultはコンソール上に出力される。
dput(df) # 出力結果をコピーしてRオブジェクトに代入すると、元と全く同一のオブジェクトが生成する。
## structure(list(face = c("ace", "two", "six"), suit = c("clubs", 
## "hearts", "diamonds"), value = c(1, 2, 3)), class = "data.frame", row.names = c(NA, 
## -3L))

1.3.10 tibble

tibbleはデータフレームの進化版。データフレームとは異なり、任意のRオブジェクト（データフレームやグラフを含む）を格納できる。 また入力データの型を自動変換しない。tidyverseではtibbleを標準の形式としている。 古い関数ではtibbleが使えないことがあるのでその場合はデータフレームを使用する。

as_tibble(iris)  # as_tibble()でdataframeからtibbleへ変換できる
## # A tibble: 150 × 5
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           <dbl>       <dbl>        <dbl>       <dbl> <fct>  
##  1          5.1         3.5          1.4         0.2 setosa 
##  2          4.9         3            1.4         0.2 setosa 
##  3          4.7         3.2          1.3         0.2 setosa 
##  4          4.6         3.1          1.5         0.2 setosa 
##  5          5           3.6          1.4         0.2 setosa 
##  6          5.4         3.9          1.7         0.4 setosa 
##  7          4.6         3.4          1.4         0.3 setosa 
##  8          5           3.4          1.5         0.2 setosa 
##  9          4.4         2.9          1.4         0.2 setosa 
## 10          4.9         3.1          1.5         0.1 setosa 
## # ℹ 140 more rows
# as.data.frame()でさらにデータフレームに戻す、データフレームは全ての行を表示してしまうので、head()を使用
head(as.data.frame(as_tibble(iris)))
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

tibbleの例、tibbleはdefaultで上から10行までしか表示しない。

library(tidyverse)
library(nycflights13)
data("flights")
flights
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
# intは整数を表す
# dblはdoubleの略で実数を表す
# chrは文字ベクトルすなわち文字列を表す
# dttmはdate-times（日付+時刻）の略
# lglはlogicalの略で、TRUEまたはFALSEからなるベクトルを示す
# fctはfactorの略で、規定の値のカテゴリ変数をRで扱う時に使用
# dateは日付を表す

tibbleの列にリストが格納されている場合、その列を特にリスト列と呼ぶ。データフレーム形式ではリストは列として扱われ、含まれるベクトルの長さが異なるとエラーになる。tibble形式ではそのような制約は無い。

data.frame(x = list(1:3, 3:5))
##   x.1.3 x.3.5
## 1     1     3
## 2     2     4
## 3     3     5
tibble(
  x = list(1:3, 3:5, 5:10),
  y = c("1, 2", "3, 4, 5", NA)
)
## # A tibble: 3 × 2
##   x         y      
##   <list>    <chr>  
## 1 <int [3]> 1, 2   
## 2 <int [3]> 3, 4, 5
## 3 <int [6]> <NA>

1.3.11 要素の選択

1.3.11.1 ベクトルの要素

角括弧[ ]内に添え字を記述して指定。添字は正の整数、負の整数、ゼロ、論理値、名前が使用できる。

vec_num <- c(1:30)
vec_num[7] # 10番目の要素、Rでは1から添え字が始まる
## [1] 7
vec_num[-3] # 3番目の要素以外
##  [1]  1  2  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## [26] 27 28 29 30
vec_num[5:13] # :で範囲を指定
## [1]  5  6  7  8  9 10 11 12 13
vec_num[-(3:10)] # 3〜10番目以外
##  [1]  1  2 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
vec_num[-3:-10] # こう書いても同じ
##  [1]  1  2 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
vec_fruits <- c(1, 2, 3)
names(vec_fruits) <- c("apple", "banana", "orange") # 各要素の名前を付ける
vec_fruits["banana"] # 名前を指定して要素を抽出
## banana 
##      2

1.3.11.2 データフレーム中の要素

ベクトルと同様、角括弧[ ]を使用。カンマで区切られた2つの添え字を記述して指定。

deck <- read.csv("deck.csv", header =T, sep = ",", stringsAsFactors = F) # csvファイルから読み込み、下記参照。
deck
##     face     suit value
## 1   king   spades    13
## 2  queen   spades    12
## 3   jack   spades    11
## 4    ten   spades    10
## 5   nine   spades     9
## 6  eight   spades     8
## 7  seven   spades     7
## 8    six   spades     6
## 9   five   spades     5
## 10  four   spades     4
## 11 three   spades     3
## 12   two   spades     2
## 13   ace   spades     1
## 14  king    clubs    13
## 15 queen    clubs    12
## 16  jack    clubs    11
## 17   ten    clubs    10
## 18  nine    clubs     9
## 19 eight    clubs     8
## 20 seven    clubs     7
## 21   six    clubs     6
## 22  five    clubs     5
## 23  four    clubs     4
## 24 three    clubs     3
## 25   two    clubs     2
## 26   ace    clubs     1
## 27  king diamonds    13
## 28 queen diamonds    12
## 29  jack diamonds    11
## 30   ten diamonds    10
## 31  nine diamonds     9
## 32 eight diamonds     8
## 33 seven diamonds     7
## 34   six diamonds     6
## 35  five diamonds     5
## 36  four diamonds     4
## 37 three diamonds     3
## 38   two diamonds     2
## 39   ace diamonds     1
## 40  king   hearts    13
## 41 queen   hearts    12
## 42  jack   hearts    11
## 43   ten   hearts    10
## 44  nine   hearts     9
## 45 eight   hearts     8
## 46 seven   hearts     7
## 47   six   hearts     6
## 48  five   hearts     5
## 49  four   hearts     4
## 50 three   hearts     3
## 51   two   hearts     2
## 52   ace   hearts     1
deck[1, 2] # 1行2列目の要素を返す
## [1] "spades"
deck[1, 1:3] # 1行1〜3列、1:3をc(1,2,3)にしても同じ
##   face   suit value
## 1 king spades    13
deck[c(2, 6:11), 2:3]  # 2行目と6〜11行目、2〜3列
##      suit value
## 2  spades    12
## 6  spades     8
## 7  spades     7
## 8  spades     6
## 9  spades     5
## 10 spades     4
## 11 spades     3
deck[-(2:40), 1:3] # 2〜40行目以外
##     face   suit value
## 1   king spades    13
## 41 queen hearts    12
## 42  jack hearts    11
## 43   ten hearts    10
## 44  nine hearts     9
## 45 eight hearts     8
## 46 seven hearts     7
## 47   six hearts     6
## 48  five hearts     5
## 49  four hearts     4
## 50 three hearts     3
## 51   two hearts     2
## 52   ace hearts     1
deck[1, ] # スペースで列全体を取り出す
##   face   suit value
## 1 king spades    13
deck[ , 2] # スペースで行全体を取り出す
##  [1] "spades"   "spades"   "spades"   "spades"   "spades"   "spades"  
##  [7] "spades"   "spades"   "spades"   "spades"   "spades"   "spades"  
## [13] "spades"   "clubs"    "clubs"    "clubs"    "clubs"    "clubs"   
## [19] "clubs"    "clubs"    "clubs"    "clubs"    "clubs"    "clubs"   
## [25] "clubs"    "clubs"    "diamonds" "diamonds" "diamonds" "diamonds"
## [31] "diamonds" "diamonds" "diamonds" "diamonds" "diamonds" "diamonds"
## [37] "diamonds" "diamonds" "diamonds" "hearts"   "hearts"   "hearts"  
## [43] "hearts"   "hearts"   "hearts"   "hearts"   "hearts"   "hearts"  
## [49] "hearts"   "hearts"   "hearts"   "hearts"
deck[ , "suit"] # オブジェクトが名前を持っているならば名前で指定可能
##  [1] "spades"   "spades"   "spades"   "spades"   "spades"   "spades"  
##  [7] "spades"   "spades"   "spades"   "spades"   "spades"   "spades"  
## [13] "spades"   "clubs"    "clubs"    "clubs"    "clubs"    "clubs"   
## [19] "clubs"    "clubs"    "clubs"    "clubs"    "clubs"    "clubs"   
## [25] "clubs"    "clubs"    "diamonds" "diamonds" "diamonds" "diamonds"
## [31] "diamonds" "diamonds" "diamonds" "diamonds" "diamonds" "diamonds"
## [37] "diamonds" "diamonds" "diamonds" "hearts"   "hearts"   "hearts"  
## [43] "hearts"   "hearts"   "hearts"   "hearts"   "hearts"   "hearts"  
## [49] "hearts"   "hearts"   "hearts"   "hearts"
deck[16, "face"]
## [1] "jack"
deck[1:10, 1, drop=F] # 1列だけ抽出する場合、データフレームがほしい場合はdrop =Fとする
##     face
## 1   king
## 2  queen
## 3   jack
## 4    ten
## 5   nine
## 6  eight
## 7  seven
## 8    six
## 9   five
## 10  four
class(deck[1:10, 1]) # 普通はベクトル型になる
## [1] "character"

ドル記号$はデータフレームやリストから名前の付いている要素をベクトルやデータフレームとして取り出すことができる。

deck$face  # $でデータフレームから列を抽出
##  [1] "king"  "queen" "jack"  "ten"   "nine"  "eight" "seven" "six"   "five" 
## [10] "four"  "three" "two"   "ace"   "king"  "queen" "jack"  "ten"   "nine" 
## [19] "eight" "seven" "six"   "five"  "four"  "three" "two"   "ace"   "king" 
## [28] "queen" "jack"  "ten"   "nine"  "eight" "seven" "six"   "five"  "four" 
## [37] "three" "two"   "ace"   "king"  "queen" "jack"  "ten"   "nine"  "eight"
## [46] "seven" "six"   "five"  "four"  "three" "two"   "ace"
typeof(deck$value)
## [1] "integer"
list1 <- list(100:130, "R", list(T, F))
names(list1) <- c("sequence", "character", "logical")  ##  namesで名前属性をつけることができる
list1$sequence  # ベクトル形式で取り出される
##  [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
typeof(list1$sequence)
## [1] "integer"
typeof(list1[1])  # 要素番号で取り出すとリスト形式になる
## [1] "list"

重要：リストの要素に名前が付いていない場合、2重の角括弧で要素を抽出すると、$記法と同じくベクトル/データフレーム形式になる
データフレーム形式になると角括弧[ ]を使用して、[行, 列]で要素が指定できる

list1[2]
## $character
## [1] "R"
typeof(list1[2])
## [1] "list"
list1[[2]]
## [1] "R"
typeof(list1[[2]])
## [1] "character"
list1$character
## [1] "R"
typeof(list1$character)
## [1] "character"
deck_list <- list(deck$face, deck$suit, deck)
deck_list[[3]][1,1,1] # データフレームなので[]で要素が指定できる
## [1] "king"
# deck_list[3][1,1,1] # deck_list[3]はリスト形式なのでエラーになる

attach()関数でデータフレームを読み込むと変数名だけで値にアクセスできるようになる。解除するにはdetach()を使う。

df <- data.frame(face=c("ace","two","six"), suit=c("spades","clubs","hearts"),  # names, class, row.names属性をもつ
                 value=c(1,2,3), stringsAsFactors = F)
df$suit
## [1] "spades" "clubs"  "hearts"
attach(df)
suit
## [1] "spades" "clubs"  "hearts"
detach(df)

attach()と似た関数にwith()がある。構文はwith(data, expr, ...)。with関数内に限ってattach(data)した後と同様にexprが評価される。

df <- data.frame(face=c("queen","ace","six"), suit=c("diamonds","spades","clubs"),  # names, class, row.names属性をもつ
                 value=c(1,2,3), stringsAsFactors = F)
with(df, paste(face, suit, sep = " of "))
## [1] "queen of diamonds" "ace of spades"     "six of clubs"

1.3.11.3 値の書き換え

1.3.11.3.1 比較演算子

式が真（TRUE）か偽（FALSE）かを判定するための演算子。式のそれぞれにTRUEかFALSEを返す。
ベクトルを比較する演算子を使った場合、Rは要素ごとに判定する。

1 < 2 # 左辺（LHS）は右辺（RHS）より小さい
## [1] TRUE
2 <= 2 # LHSはRHSと等しいか小さい
## [1] TRUE
5 > 6  # LHSはRHSより大きい
## [1] FALSE
4 >= 5 # LHSはRHSと等しいか大きい
## [1] FALSE
(2 + 3) == (4 + 1)
## [1] TRUE
((2 * 3) + 1) != (2 * (3 + 1))
## [1] TRUE
1 > 0:2
## [1]  TRUE FALSE FALSE
c(1, 2, 3) == c(3, 2, 1)
## [1] FALSE  TRUE FALSE

マッチ演算子 %in% は左側の個々の値が右側のベクトルに含まれているかどうかをテストする。
従って左側の要素の個数だけTRUEまたはFALSEを返す。

1 %in% c(3, 4, 5)
## [1] FALSE
c(1, 2) %in% c(3, 4, 5)
## [1] FALSE FALSE
c(1, 2, 3, 4) %in% c(3, 4, 5)
## [1] FALSE FALSE  TRUE  TRUE
1.3.11.3.2 ブール演算子

ブール演算子は複数の論理テストから1つのTRUEかFALSEを導出する。Rは個々の論理テストを実行してから、ブール演算子を使用して判定する。ベクトルとともに使うと、ブール演算子は要素単位で評価する。

a <- c(1, 2, 3)
b <- c(1, 2, 3)
c <- c(1, 2, 4)
a == b
## [1] TRUE TRUE TRUE
b == c
## [1]  TRUE  TRUE FALSE
!(b == c) # 論理テストの反転
## [1] FALSE FALSE  TRUE
a == b & b == c  # cond1 & cond2, cond1とcond2が両方ともTRUEかどうか
## [1]  TRUE  TRUE FALSE
a == b | b == c  # cond1 & cond2, cond1とcond2のどちらかがTRUEかどうか
## [1] TRUE TRUE TRUE
any(a == b, b == c)  # any(cond1, cond2, ...) condの中で1つでもTRUEになっているものがあるか
## [1] TRUE
all(a == b, b == c)  # all(cond1, cond2,...) 全てのcondがTRUEかどうか
## [1] FALSE

ダブル演算子&&や||はそれぞれ&や|と同じように動作するが、ベクトル対応になっておらずスカラーに対して適用され、結果はTRUEまたはFALSE。

# a == b && b == c
# a == b || b == c
1.3.11.3.3 論理テストによる要素の抽出と書き換え

ベクトルの要素はTRUE/FALSEの配列で指定できる。

vec_num <- sample(1:6, 6, replace = FALSE)  # 1:6をランダムに6個並べた配列を作る
vec_num
## [1] 2 1 3 6 4 5
vec_num[c(FALSE, TRUE, FALSE, TRUE, TRUE, FALSE)] # 2, 4, 5番目の要素を抽出
## [1] 1 6 4

従って論理テストの結果を利用して、要素を抽出することができる。sum()関数内でTRUEとFALSEはそれぞれ1と0に置換して計算されるので、条件に当てはまる要素の個数を数えることができる。また代入演算子を使用して値を置換することも可能。

vec_num2 <- sample(1:6, 20, replace = TRUE)  # 1:6をランダムに20個並べた配列を作る
vec_num2
##  [1] 6 3 3 1 2 3 4 3 2 2 6 5 5 4 3 6 5 6 6 2
(vec_num2 %% 2) == 0  # 要素が2で割り切れる（=偶数）かどうかを判定
##  [1]  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE
## [13] FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE
vec_num2[(vec_num2 %% 2) == 0] # 偶数だけを抽出
##  [1] 6 2 4 2 2 6 4 6 6 6 2
vec_num2 %in% c(2, 5) # 要素が2または5かどうかを判定
##  [1] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE
## [13]  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE
sum(vec_num2 %in% c(2, 5)) # 2または5の個数をカウント
## [1] 7
vec_num2[vec_num2 %in% c(2, 5)] # 2または5を抽出
## [1] 2 2 2 5 5 5 2
vec_num2[vec_num2 %in% c(2, 5)] <- 0 # 2または5を0に置換
vec_num2
##  [1] 6 3 3 1 0 3 4 3 0 0 6 0 0 4 3 6 0 6 6 0

データフレームの場合も同様に論理テストが利用できる。

deck <- read.csv("deck.csv", header =T, sep = ",", stringsAsFactors = F) # csvファイルから読み込み
head(deck)
##    face   suit value
## 1  king spades    13
## 2 queen spades    12
## 3  jack spades    11
## 4   ten spades    10
## 5  nine spades     9
## 6 eight spades     8
sum(deck$face == "ace")  # aceの個数を数える
## [1] 4
deck$value[deck$face == "queen"]  # queenのポイントを表示
## [1] 12 12 12 12
deck[deck$suit == "spades", ]  # spadesのカードを全て表示
##     face   suit value
## 1   king spades    13
## 2  queen spades    12
## 3   jack spades    11
## 4    ten spades    10
## 5   nine spades     9
## 6  eight spades     8
## 7  seven spades     7
## 8    six spades     6
## 9   five spades     5
## 10  four spades     4
## 11 three spades     3
## 12   two spades     2
## 13   ace spades     1
deck$suit == "hearts" & deck$face == "king"  # ハートのキングを判定する論理テスト
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE
kingOfHearts <- deck$suit == "hearts" & deck$face == "king" # テストの結果をオブジェクトに格納
deck[kingOfHearts, ]  # ハートのキングのカードを表示
##    face   suit value
## 40 king hearts    13
facecards <- deck$face %in% c("king", "queen", "jack")  # フェイスカード（キング、クイーン、ジャック）を判定
deck[facecards, ]
##     face     suit value
## 1   king   spades    13
## 2  queen   spades    12
## 3   jack   spades    11
## 14  king    clubs    13
## 15 queen    clubs    12
## 16  jack    clubs    11
## 27  king diamonds    13
## 28 queen diamonds    12
## 29  jack diamonds    11
## 40  king   hearts    13
## 41 queen   hearts    12
## 42  jack   hearts    11
deck$value[facecards]
##  [1] 13 12 11 13 12 11 13 12 11 13 12 11

1.4 関数の作成

コードの塊を3回以上コピー&ペーストする場合、関数化を考えるべき。
関数を書くことのメリット
* コードが理解しやすい直感的な名前を関数につけることができる。 * 要求が変わっても、複数箇所を変更せずに、1ヶ所のコードを変更するだけで済む
* コピー&ペースト時の間違いをなくすことができる

df <- tibble::tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

df$a <- (df$a - min(df$a, na.rm = TRUE)) / 
  (max(df$a, na.rm = TRUE) - min(df$a, na.rm = TRUE))
df$b <- (df$b - min(df$b, na.rm = TRUE)) / 
  (max(df$b, na.rm = TRUE) - min(df$a, na.rm = TRUE))  ## 間違いあり
df$c <- (df$c - min(df$c, na.rm = TRUE)) / 
  (max(df$c, na.rm = TRUE) - min(df$c, na.rm = TRUE))
df$d <- (df$d - min(df$d, na.rm = TRUE)) / 
  (max(df$d, na.rm = TRUE) - min(df$d, na.rm = TRUE))

(df$a - min(df$a, na.rm = TRUE)) /
  (max(df$a, na.rm = TRUE) - min(df$a, na.rm = TRUE))
##  [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
##  [7] 1.00000000 0.64768651 0.79237052 0.18434257
x <- df$a
(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
##  [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
##  [7] 1.00000000 0.64768651 0.79237052 0.18434257
#>  [1] 0.2892677 0.7509271 0.0000000 0.6781686 0.8530656 1.0000000 0.1716402
#>  [8] 0.6107464 0.6116181 0.6008793
rng <- range(x, na.rm = TRUE)  ## rangeは最小値と最大値を含んだベクトルを返す
(x - rng[1]) / (rng[2] - rng[1])  ## rng[1]は最小値、rng[2]は最大値
##  [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
##  [7] 1.00000000 0.64768651 0.79237052 0.18434257

関数化のステップ:1. 関数に名前をつける, 2. 引数を関数の内側に書く, 3. 関数直後の{}にコードを書く。

rescale01 <- function(x) {  
rng <- range(x, na.rm = TRUE)  
(x - rng[1]) / (rng[2] - rng[1])
}
rescale01(c(0, 5, 10)) # 他の値で上手くいくか確認
## [1] 0.0 0.5 1.0

1.4.1 if文による条件実行

ベクトルの各要素に名前があるかどうかを示す論理ベクトルを返す関数。戻り値を指定しない場合、最後に計算した値を返す。

has_name <- function(x){
  nms <- names(x)
  if (is.null(nms)){
    rep(FALSE, length(x))
  } else{
      !is.na(nms) & nms != ""
    }
}

|| (or)や&& (and)を使用して複数の論理式を結合できる。
||は最初にTRUEと評価された要素でTRUEを返し、残りの要素は評価しない。
&&ではFALSEと評価した最初の要素で、FALSEを返す。
if文では|や&を使用してはいけない、これらは複数値に対するベクトル演算を行う（なのでfilter()で使用する）。
論理ベクトルはany()やall()を使用して、単一値にまとめることができる。
等しいかどうかのテスト == もベクトル演算なので注意する。

a=c(1,1,1,1,1); b=c(1,1,0,0,1)
a == b
## [1]  TRUE  TRUE FALSE FALSE  TRUE
identical(a, b)  ## identical は非常に厳格で常に単一のTかFを返し、型強制はしない
## [1] FALSE

1.4.2 複合条件

複数のif文をつなげることができる：
if (this) {
# do that
} else if (that) {
# do something else
} else {
}
switch関数を使用すると、位置や名前に基づいて選択したコードを評価できる。

op <- "1"
x <- 10
y <- 20
ans <- switch(op,
              "1" = x + y,
              "2" = x - y,
              "3" = x * y,
              "4" = x / y,
              stop("Only can use 1, 2, 3, and 4"))
ans
## [1] 30

1.4.3 関数の引数

一つは計算するデータ、もう一つは計算の詳細を制御。例えばlog()ではデータがx、詳細が対数の底base。一般にデータ引数が最初に来る、詳細引数は最後に来て、通常はデフォルト値が存在。

# 正規近似を使って、平均値の信頼区間を計算
mean_ci <- function(x, conf = 0.95) { ## 0.95がdefaultの値
  se <- sd(x) / sqrt(length(x))
  alpha <- 1 - conf
  mean(x) + se * qnorm(c(alpha / 2, 1 - alpha / 2))
}
x <- runif(100)
mean_ci(x)
## [1] 0.4164526 0.5251495
#> [1] 0.4976111 0.6099594
mean_ci(x, conf = 0.99)
## [1] 0.3993751 0.5422271

1.4.4 値のチェック

重み付き要約統計量を計算する関数

wt_mean <- function(x, w) {
  sum(x * w) / sum(w)
}
wt_var <- function(x, w) {
  mu <- wt_mean(x, w)
  sum(w * (x - mu) ^ 2) / sum(w)
}
wt_sd <- function(x, w) {
  sqrt(wt_var(x, w))
}

wt_mean(1:6, 1:5)  ## xとwのベクトルの長さが異なるとエラーになる
## Warning in x * w:
## 長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません
## [1] 4.066667

重要な前提条件はチェックして、正しくないとstopと共にエラーを返すようにする

wt_mean <- function(x, w) {
  if (length(x) != length(w)) {
    stop("`x` and `w` must be the same length", call. = FALSE)  ## xとwの長さが異なる場合、stopとともにエラーメッセージを返す
  }
  sum(w * x) / sum(w)
}

stopifnot()は各引数がTRUEかどうかを調べ、そうでないと汎用的なエラーメッセージを返す

wt_mean <- function(x, w, na.rm = FALSE) {
  stopifnot(is.logical(na.rm), length(na.rm) == 1)
  stopifnot(length(x) == length(w))
  
  if (na.rm) {
    miss <- is.na(x) | is.na(w)
    x <- x[!miss]
    w <- w[!miss]
  }
  sum(w * x) / sum(w)
}
# wt_mean(1:6, 6:1, na.rm = "foo")

1.4.5 3ドット

Rの多くの関数では入力する引数の数に制限は無い、例えばsumやstr_cなど。

sum(1,2,3,4,5,6,7,8,9,10)
## [1] 55
stringr::str_c("a","b","c","d","e","f","g")
## [1] "abcdefg"

特別な引数「…」は入力引数がいくつでもマッチでき、他の関数に送り出すことができる。

commas <- function(...) stringr::str_c(..., collapse = ", ")
commas(letters[1:10])
## [1] "a, b, c, d, e, f, g, h, i, j"
rule <- function(..., pad = "-") {
  title <- paste0(...)
  width <- getOption("width") - nchar(title) - 5   ## getOption("width")はコンソール画面の幅
  cat(title, " ", stringr::str_dup(pad, width), "\n", sep = "")
}
rule("Important output")
## Important output -----------------------------------------------------------

1.4.6 戻り値

関数が返す値は通常は評価の最後の文。return()を使って早く返すように選択できる。考慮するべきは1. 値を早く返すと関数が読みやすくなるか 2. パイプを使って関数が書けるか。 明示的return文にするべきは主に次の3つのケース： 1. 入力が空の場合、2. 複雑なブロックが一つと簡単なブロックが一つある場合、3. パイプにできる関数を書く場合

1.5 イテレーション

重複をなくすツールの1つが関数化で、繰り返し使われるコードパターンを独立させて、たやすく再利用や更新できるようにする。もう１つのツールがイテレーション（iteration）で、複数の入力（異なる列やデーセット）に対して同じことを行う。イテレーションには命令型（for loop, while loop）と関数型がある。

1.5.1 forループ

forループの3つの要素:
出力: ループを開始する前に、出力用のスペースを確保しておく必要がある。指定長の空ベクトルを作るにはvector()を使う。
シーケンス: 何でループするのかを決める。seq_along()は1からベクトルの長さ（行数）までの数列を生成、ベクトルの長さが0の時にもうまく処理してくれるので1:length()よりも安全。
本体: 実際の仕事をするコード。

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)
# for loopで計算
output <- vector("double", ncol(df))  ## 1. 出力、結果を出力するための指定長の空ベクトルを作っておく
for (i in seq_along(df)){  ## 2. シーケンス
  output[[i]] <- median(df[[i]])  ## 3. 本体
}
output
## [1] -0.0962729  0.3255362 -0.4405117 -0.6236106

[[]]は単一要素を扱うことが明示されるので、for loopの中では[[]]を使用する（例えアトミックベクトルでも[[]]の使用を推奨）。

1.5.2 forループのバリエーション

forループには以下の4つのバリエーションがある：

  • 新たなオブジェクトを作らず、既存オブジェクトを変更
  • 添字ではなく、名前や値についてループする
  • 長さが不明な出力を扱う
  • 長さの不明なシーケンスを扱う 

1.5.2.1 既存オブジェクトの変更

出力は既にある、すなわち入力と同じ。データフレームを列のリストと考え、列をiteration。

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)
# 関数で解く場合
rescale01 <- function(x) {
  rng <- range(x, na.rm = TRUE)
  (x - rng[1]) / (rng[2] - rng[1])
}

# for loopで解く場合
for (i in seq_along(df)){ 
  df[[i]] <- rescale01(df[[i]])  ## 本体：rescale01()する
}

1.5.2.2 ループパターン

ベクトルでループするには基本的に3つの方法がある：

  • 添字でループ：最も一般的な方法、(i in seq_along(xs))で数値の添字を使って値xs[[i]]を抽出。names(xs)[[i]]で名前も取得できる。
  • 要素を使ってループ： for (x in xs)、出力を効率的に格納するのは難しい。グラフを書いたりするなど、副作用だけを使う際に有効。
  • 名前を使ってループ：(nm in names(xs))、xs[[nm]]で値にアクセスできる名前を指定。名前付きの出力を作るなら、結果ベクトルに名前を付けるのを忘れないようにする。
results <- vector("list", length(x))
names(results) <- names(x)

1.5.2.3 出力長不明

出力の長さがどれだけになるかわからない場合もある。例：長さがランダムなベクトルをシミュレーションしたい。

means <- c(0, 1, 2)
output <- double()
for (i in seq_along(means)) {
  n <- sample(100, 1)  ## 1〜100の間の整数を一つランダムに選ぶ
#ベクトルoutputの最後にrnorm(..)で選んだn個の数字を追加して、ベクトルを長くしていく
  output <- c(output, rnorm(n, means[[i]]))  
}
str(output)
##  num [1:152] -2.255 -0.864 0.621 0.143 0.951 ...

これでは毎回Rが前のデータを全てコピーしないといけないので、あまり効率的ではない。 より良い解法では、結果をリストに格納し、ループが終わってから一つのベクトルにまとめる。

out <- vector("list", length(means))  ## 長さが3のリストを作成できる
for (i in seq_along(means)) {
  n <- sample(100, 1) 
  out[[i]] <- rnorm(n, means[[i]])  ## リストのi番目の要素に格納
}
str(out)
## List of 3
##  $ : num [1:10] -0.594 1.064 1.116 -0.397 1.727 ...
##  $ : num [1:93] 1.5048 1.8186 -0.0516 -0.3624 0.7088 ...
##  $ : num [1:92] 3.31 2.88 3.91 2.02 1.27 ...
str(unlist(out))  ## ベクトルのリストを1つのベクトルにする、purrr::flatten_dbl()を使う方がベター
##  num [1:195] -0.594 1.064 1.116 -0.397 1.727 ...

他の例:

  • 長い文字列を生成する場合、出力をベクトルに保存しておき、paste(output, collapse="")でベクトルの要素を連結して文字列にする。
  • 大きなデータフレームを生成する場合、繰り返しのたびにrbind()で連結するのではなく、出力をリストに保存し、dplyr::bind_rows(output)を使って連結して１つのデータフレームにまとめる。

1.5.2.4 シーケンス長不明

何回繰り返すか不明の場合はwhileを使用。例えば表が3回連続出るまでに何回硬貨投げを試行するのかを数えるwhileループ。

flip <- function() sample(c("T", "H"), 1)
flips <- 0
nheads <- 0

while (nheads < 3) {
  if (flip() == "H") {
    nheads <- nheads + 1
  } else {
    nheads <- 0
  }
  flips <- flips + 1
}
flips
## [1] 42

1.5.2.5 forループと関数型

Rは関数型プログラミング言語であり、他の言語ほどforループが重要とは考えない。forループは関数でラップ（wrap）できる。

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)
# for loopで各列の平均を計算
output <- vector("double", length(df))
for (i in seq_along(df)) {
  output[[i]] <- mean(df[[i]])
}
output
## [1] -0.2607816 -0.8224295  0.1765180 -0.3967573

平均値の計算を関数化。中央値と標準偏差の計算も同様に関数化。

col_mean <- function(df) {   ## 平均の計算を関数化
  output <- vector("double", length(df))
  for (i in seq_along(df)) {
    output[i] <- mean(df[[i]])
  }
  output           # 返り値
  # return(output) # こう書いても同じ
}
col_mean(df)
## [1] -0.2607816 -0.8224295  0.1765180 -0.3967573
col_median <- function(df) { ## 中央値の計算を関数化
  output <- vector("double", length(df))
  for (i in seq_along(df)) {
    output[i] <- median(df[[i]])
  }
  output
}
col_median(df)
## [1] -0.2797751 -0.5692534  0.2601980 -0.4112665
col_sd <- function(df) {  #### 標準偏差の計算を関数化
  output <- vector("double", length(df))
  for (i in seq_along(df)) {
    output[i] <- sd(df[[i]])
  }
  output
}
col_sd(df)
## [1] 0.4467616 1.0601581 0.4877669 1.0318189

平均値、中央値、標準偏差の関数を引数として渡せば、より一般的した関数を定義できる。関数を引数にする場合は()を付けない。

# 関数を別の関数に引き渡す
col_summary <- function(df, fun) {  ## 引数funを追加、mean, median, sd等を入力すれば良い
  out <- vector("double", length(df))
  for (i in seq_along(df)) {
    out[i] <- fun(df[[i]])
  }
  out
}
col_summary(df, median)
## [1] -0.2797751 -0.5692534  0.2601980 -0.4112665
col_summary(df, mean)
## [1] -0.2607816 -0.8224295  0.1765180 -0.3967573

1.5.3 マップ関数

purrrパッケージのマップ関数とは第1引数のベクトルの各要素について何かを行い、結果を保存するための関数。どの関数もベクトルを入力し、各要素に関数を適用し、入力と同じ長さ（同じ名前）の新たなベクトルを返す。第2引数.fは式、文字ベクトル、または整数ベクトルでも良い。3ドット「…」で.fが呼ばれる際の引数を渡す。文字や整数の場合はその位置の要素が抽出される（後述）：

  • map(): makes a list., lapplyと同様
  • map_lgl(): makes a logical vector.
  • map_int(): makes an integer vector.
  • map_dbl(): makes a double vector.
  • map_chr(): makes a character vector.
map(df, mean) # リストを返す、各列の名前も保持される
## $a
## [1] -0.2607816
## 
## $b
## [1] -0.8224295
## 
## $c
## [1] 0.176518
## 
## $d
## [1] -0.3967573
map_dbl(df, mean)  # ベクトルを返す
##          a          b          c          d 
## -0.2607816 -0.8224295  0.1765180 -0.3967573
map_chr(df, mean)  # 要素が文字型のベクトルを返す
## Warning: Automatic coercion from double to character was deprecated in purrr 1.0.0.
## ℹ Please use an explicit call to `as.character()` within `map_chr()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
##           a           b           c           d 
## "-0.260782" "-0.822429"  "0.176518" "-0.396757"
map_dbl(df, median)
##          a          b          c          d 
## -0.2797751 -0.5692534  0.2601980 -0.4112665
map_dbl(df, sd)
##         a         b         c         d 
## 0.4467616 1.0601581 0.4877669 1.0318189
df %>% map_dbl(max) # パイプ演算子も使える
##         a         b         c         d 
## 0.5653385 0.5763737 0.7829316 1.1426596

1.5.3.1 ショートカット

mtcarsデータセットをシリンダーの値で3つに分け、同じ線形モデルを適用したい。

models <- mtcars %>% 
  split(.$cyl) %>% 
  map(function(df) lm(mpg ~ wt, data = df))  ## 匿名関数を定義して使用

匿名関数の代わりに、one-sided fomulaオブジェクト（~）を使用して次のように書ける。ここでピリオドは現在のリスト要素を参照する代名詞のようなもの。

models <- mtcars %>% 
  split(.$cyl) %>% 
  map(~ lm(mpg ~ wt, data = .))

R^2のような要約統計量を取得するにはsummary()を実行してから、r.squaredという成分を抽出する。

models %>%
  map(summary) %>%   ## summaryを実行してから、リストに格納されているR2値成分を抽出する
  map_dbl(~ .$r.squared)
##         4         6         8 
## 0.5086326 0.4645102 0.4229655
# 文字でも抽出できる、こちらの方が一般的
models %>% 
  map(summary) %>% 
  map_dbl("r.squared")
##         4         6         8 
## 0.5086326 0.4645102 0.4229655

整数を使ってその位置の要素を抽出することもできる。

x <- list(list(1, 2, 3), list(4, 5, 6), list(7, 8, 9))
x %>% map_dbl(2)
## [1] 2 5 8

1.5.3.2 apply系関数とpurrr::map関数の類似性

  • base系のlapply()は基本的にmap()と同じ。map()の方がpurrrの他の関数との一貫性があり、.fのショートカットを使う。
  • sapply()はlapply()のwrapperで、出力を自動的に簡素化。どんな出力になるかわからないので自作関数に使うには問題がある。
  • vapply()は型を定義する引数を追加するので、sapply()の安全版。入力が面倒になる。vapply(df, is.numeric, logical(1))はmap_lgl(df, is.numeric)と等価。map関数はベクトルしか作成できないのに対し、vapply()は行列も作成できる。

map関数の練習問題：

# a. mtcarsの各列の平均を計算する
data("mtcars")
mtcars %>% map_dbl(mean)
##        mpg        cyl       disp         hp       drat         wt       qsec 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
##         vs         am       gear       carb 
##   0.437500   0.406250   3.687500   2.812500
# b. nycflights13::flightsの各列の型を決定する
data("flights")
flights %>% map_chr(typeof)
##           year          month            day       dep_time sched_dep_time 
##      "integer"      "integer"      "integer"      "integer"      "integer" 
##      dep_delay       arr_time sched_arr_time      arr_delay        carrier 
##       "double"      "integer"      "integer"       "double"    "character" 
##         flight        tailnum         origin           dest       air_time 
##      "integer"    "character"    "character"    "character"       "double" 
##       distance           hour         minute      time_hour 
##       "double"       "double"       "double"       "double"
# c. irisの各列の一意な値の個数を計算する
data("iris")
iris %>% map(unique) %>% map(length)
## $Sepal.Length
## [1] 35
## 
## $Sepal.Width
## [1] 23
## 
## $Petal.Length
## [1] 43
## 
## $Petal.Width
## [1] 22
## 
## $Species
## [1] 3
# 以下のようにより簡潔に書ける
iris %>% map_int(n_distinct)
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
##           35           23           43           22            3
# d. mu = -10, 0, 10, 100のそれぞれについて10個の正規乱数を生成する
map(c(-10,0,10,100), rnorm, n = 10)  ## rnorm()の引数は後ろにつける
## [[1]]
##  [1]  -9.693572  -9.105398  -9.972139 -11.043093  -9.622109  -9.986116
##  [7] -10.782066  -9.902077 -11.154659  -8.897776
## 
## [[2]]
##  [1] -1.70493586 -0.11161431 -0.29043668  0.44228369 -1.74003872 -1.31723386
##  [7]  0.03207093  0.64853885 -0.36374951  1.66519574
## 
## [[3]]
##  [1] 10.192985  9.696792  8.768456  9.734266  8.451277  9.981873 10.212536
##  [8]  9.970477 11.647382 10.735504
## 
## [[4]]
##  [1] 100.45297 102.56478  99.23805  98.68471 100.63526  99.39195  98.66007
##  [8] 100.41030 100.65494  99.92367
# 匿名関数を使用する場合
map(c(-10,0,10,100), function(mu) rnorm(n = 10, mean = mu)) 
## [[1]]
##  [1] -10.692553 -11.468975  -9.678817  -8.837789 -10.170395 -10.416223
##  [7]  -8.806419  -9.002832  -9.165819 -11.525696
## 
## [[2]]
##  [1] -1.0843696 -0.4529417 -1.1100316  0.7267748  1.4531480 -1.2733976
##  [7] -0.7647024 -0.2249176  1.8354420  0.4604920
## 
## [[3]]
##  [1]  9.633834  8.752625  9.148156 10.634914 11.378216 10.467256  9.944206
##  [8]  8.456584 10.635121 10.217572
## 
## [[4]]
##  [1]  99.19464 100.77860 100.11078  97.18729  99.78555  99.02519 101.85769
##  [8] 100.57560  99.61918  99.98137
# one-sided formulaオブジェクトを使用する場合
map(c(-10, 0, 10, 100), ~ rnorm(n = 10, mean = .))
## [[1]]
##  [1]  -9.182589  -9.737206  -9.824892  -8.893812  -9.167568  -9.988400
##  [7] -10.590297  -9.541603 -10.193090 -12.055588
## 
## [[2]]
##  [1] -0.03976529 -0.40782951  1.04600315  0.72389716 -0.15182399  0.62638232
##  [7] -0.58032425  1.00933751 -0.55416588 -0.41401237
## 
## [[3]]
##  [1] 10.349641  9.652933 10.646024 10.034130 10.263762  9.814422 11.122525
##  [8]  9.144291  9.787309  8.294519
## 
## [[4]]
##  [1] 100.34735  99.94242 100.89763  99.04754  99.88458 101.51632 101.69118
##  [8] 100.81106 100.69948 101.17728
# 3. map(1:5, runif)は何をするか、それはなぜか
map(1:5, runif)
## [[1]]
## [1] 0.6899524
## 
## [[2]]
## [1] 0.1154215 0.5156386
## 
## [[3]]
## [1] 0.9582620 0.5123572 0.5954985
## 
## [[4]]
## [1] 0.31071324 0.99986653 0.72981389 0.02290406
## 
## [[5]]
## [1] 0.1309487 0.1608921 0.6092839 0.3920713 0.1263547
# 以下と等価。
map(seq(1:5), runif)
## [[1]]
## [1] 0.1762757
## 
## [[2]]
## [1] 0.3664654 0.7194777
## 
## [[3]]
## [1] 0.8599465 0.1532712 0.5284390
## 
## [[4]]
## [1] 0.20072910 0.18444203 0.06565291 0.73609323
## 
## [[5]]
## [1] 0.7399966 0.4658728 0.7744339 0.3840219 0.5345716
# 4. map(-2:2, rnorm, n = 5)は何をするか、map_dbl(-2:2, rnorm, n = 5)はどうなるか
map(-2:2, rnorm, n = 5)
## [[1]]
## [1] -2.0593755 -1.6782743 -0.9160276 -2.3464813  0.2129120
## 
## [[2]]
## [1] -0.5782569 -0.3170911 -2.6440907 -1.6742033 -1.3217819
## 
## [[3]]
## [1] 0.85090377 0.65089650 0.09451189 0.10839872 1.13291065
## 
## [[4]]
## [1]  1.211274 -0.739638  2.442857  2.185321  2.414876
## 
## [[5]]
## [1] 2.488894 1.520707 3.071209 2.076521 1.454397
# map_dbl(-2:2, rnorm, n = 5) # エラーになる

1.5.3.3 失敗の処理

関数safely()は関数を引数に取り、修正版を返す。この場合、修正された関数はエラーを投げない。代わりに次の2つのオブジェクトをリストにまとめて返す：
result: 元の結果。エラーがあればNULL。
error: エラーオブジェクト。処理が成功すればNULL。
関数が成功すれば、要素resultには結果が入り、要素errorはNULLになる。

safe_log <- safely(log)  ## safelyはresultとerrorを1つのリストにまとめて返す
str(safe_log(10))
## List of 2
##  $ result: num 2.3
##  $ error : NULL
str(safe_log("a"))
## List of 2
##  $ result: NULL
##  $ error :List of 2
##   ..$ message: chr " 数学関数に数値でない引数が渡されました "
##   ..$ call   : language .Primitive("log")(x, base)
##   ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"

safely()はmapで使うことができる。list(result, error)を要素とするリストが得られる。

x <- list(1, 10, "a")
y <- x %>% map(safely(log))
str(y)
## List of 3
##  $ :List of 2
##   ..$ result: num 0
##   ..$ error : NULL
##  $ :List of 2
##   ..$ result: num 2.3
##   ..$ error : NULL
##  $ :List of 2
##   ..$ result: NULL
##   ..$ error :List of 2
##   .. ..$ message: chr " 数学関数に数値でない引数が渡されました "
##   .. ..$ call   : language .Primitive("log")(x, base)
##   .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"

purrr::list_transpose()を得られたリストに対して使うと、resultとerror、それぞれのリスト（list(result), list(error)）に変換される。list_transpose()は行列転置関数t()のlist版。 例えば、pair of lists <=> list of pairs。 data.frameをlistとして渡すとlist of rowsが得られる。

y <- y %>% list_transpose()
str(y)
## List of 2
##  $ result:List of 3
##   ..$ : num 0
##   ..$ : num 2.3
##   ..$ : NULL
##  $ error :List of 3
##   ..$ : NULL
##   ..$ : NULL
##   ..$ :List of 2
##   .. ..$ message: chr " 数学関数に数値でない引数が渡されました "
##   .. ..$ call   : language .Primitive("log")(x, base)
##   .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"

possibly()はsafely()同様、常に成功する、エラーの際default値を返す。quietly()はエラーを捕捉せずに出力、メッセージ、警告を補足する。

x <- list(1, 10, "a")
x %>% map_dbl(possibly(log, NA_real_))
## [1] 0.000000 2.302585       NA
x <- list(1, -1)
x %>% map(quietly(log)) %>% str() 
## List of 2
##  $ :List of 4
##   ..$ result  : num 0
##   ..$ output  : chr ""
##   ..$ warnings: chr(0) 
##   ..$ messages: chr(0) 
##  $ :List of 4
##   ..$ result  : num NaN
##   ..$ output  : chr ""
##   ..$ warnings: chr " 計算結果が NaN になりました "
##   ..$ messages: chr(0)

1.5.3.4 複数引数へのマップ

単一の入力ではなく、複数の関連する入力について並列にイテレーションする場合、関数map2()とpmap()を使う。以下は平均値の異なる正規乱数を複数シミュレーションする例。

# 平均値の異なる正規乱数を複数シミュレーションする
mu <- list(5, 10, -3) # 試す平均値のリスト
# mu <- c(5, 10, -3)  ## ベクトルにしても結果は同じ
mu %>% 
  map(rnorm, n = 5) %>% 
  str()
## List of 3
##  $ : num [1:5] 3.8 6.25 4.66 5.32 2.8
##  $ : num [1:5] 9.38 10.44 10.27 9.18 11.62
##  $ : num [1:5] -3.13 -4.09 -1.75 -3.11 -1.48

平均値に加えて標準偏差sigmaも変化させる場合。平均と標準偏差のベクトルに添え字を当ててiterationするやり方。

sigma <- list(1, 5, 10)   
# sigma <- ｃ(1, 5, 10)   
seq_along(mu) %>% 
  map(~ rnorm(n = 5, mean = mu[[.]], sd = sigma[[.]])) %>% # n =, mean = , sd =, は省略可
  str()
## List of 3
##  $ : num [1:5] 5.37 2.84 5.35 5.93 5.74
##  $ : num [1:5] 18.32 12.18 12.21 9.17 7.93
##  $ : num [1:5] -4.833 -0.727 -5.323 0.499 1.683

引数を2つ取れる、map2()で2つのベクトルを並列に処理。コードの意図がわかりやすくなる。変化する引数は関数の前に、変化しない引数は関数の後に書くことに注意。

map2(mu, sigma, rnorm, n = 5) %>% 
  str()
## List of 3
##  $ : num [1:5] 6.33 5.4 3.72 6.22 6.15
##  $ : num [1:5] 13.46 14.33 8.86 11.68 8.82
##  $ : num [1:5] -1.06 -4.58 6.67 -9.31 6.25

map3(), map4(),..等も考えられるが、purrrには代わりに引数のリストを引数として取るpmap()が用意されている。以下は平均値、標準偏差、サンプルサイズを変化させる場合。

mu <- list(5, 10, -3) # 試す平均値のリスト
sigma <- list(1, 5, 10)  
n <- list(1, 3, 5)
args1 <- list(n, mu, sigma)
args1 %>%
  pmap(rnorm) %>% 
  str()
## List of 3
##  $ : num 4.48
##  $ : num [1:3] 7.41 10.66 7.51
##  $ : num [1:5] -19.399 7.398 4.153 0.328 -12.576

リストの要素に名前を付けないと、pmap()は位置を代わりに使うので、引数に名前を付けた方が良い。

args2 <- list(mean = mu, sd = sigma, n = n)  ## 引数に名前を付けた方がわかりやすい
args2 %>% 
  pmap(rnorm) %>% 
  str()
## List of 3
##  $ : num 4.72
##  $ : num [1:3] 8.231 4.867 -0.321
##  $ : num [1:5] -2.91 -11.074 -11.724 -10.815 -0.879

引数は全て同じ長さなので、データフレームに格納しておくとより便利でわかりやすい。

params <- tribble(
  ~mean, ~sd, ~n,
  5,     1,  1,
  10,     5,  3,
  -3,    10,  5
)
params %>% 
  pmap(rnorm)
## [[1]]
## [1] 3.769903
## 
## [[2]]
## [1] -0.7947732  9.6287552  6.5797799
## 
## [[3]]
## [1]  10.3802031  -4.3383396  -0.4663692 -10.7885070   8.9599208

1.5.3.5 様々な関数を呼び出す

もう一段階複雑にして、関数への引数だけでなく、関数も変化させるにはinvoke_map()を使用する。第1引数は関数のリスト、第2引数は関数ごとに異なる引数を指定するリストのリスト、第3引数以降は全関数に渡される。

f <- c("runif", "rnorm", "rpois")
param <- list(
  list(min = -1, max = 1), 
  list(sd = 5), 
  list(lambda = 10)
)
invoke_map(f, param, n = 5) %>% str()
## Warning: `invoke_map()` was deprecated in purrr 1.0.0.
## ℹ Please use map() + exec() instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## List of 3
##  $ : num [1:5] 0.964 -0.384 0.238 0.489 0.565
##  $ : num [1:5] -9.06 2.68 -7.17 -4.65 3.49
##  $ : int [1:5] 8 9 9 10 15

パラメータと計算結果を1つのデータフレームにまとめることもできる。

sim <- tribble(   ## データフレームに関数と引数のリストを格納
  ~f,      ~params,
  "runif", list(min = -1, max = 1),
  "rnorm", list(sd = 5),
  "rpois", list(lambda = 10)
)
sim %>% 
  mutate(sim = invoke_map(f, params, n = 10))  ## シミュレーション結果を列simに格納
## # A tibble: 3 × 3
##   f     params           sim       
##   <chr> <list>           <list>    
## 1 runif <named list [2]> <dbl [10]>
## 2 rnorm <named list [1]> <dbl [10]>
## 3 rpois <named list [1]> <int [10]>

1.5.3.6 ウォーク

ウォークは関数をその戻り値ではなく、副作用のために呼び出したい際にmapの代わりに使う。関数を適用しつつ第一引数を暗黙的に返す。グラフのリストとファイル名のベクトルがあると、pwalk()を使って、各ファイルをディスクの対応する場所に保存できる。

plots <- mtcars %>% 
  split(.$cyl) %>% 
  map(~ ggplot(., aes(mpg, wt)) + geom_point())
paths <- stringr::str_c(names(plots), ".pdf")

pwalk(list(paths, plots), ggsave, path = tempdir())
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image

1.5.3.7 forループの他のパターン

purrrにはforループの他のパターンを抽象化する関数も用意されている。
TRUEかFALSEを返す述語関数とともに機能する関数、keep()とdiscard()は述語がTRUEまたはFALSEとなる要素をそれぞれ返す。

iris %>% 
  keep(is.factor) %>% 
  str()
## 'data.frame':    150 obs. of  1 variable:
##  $ Species: Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
iris %>% 
  discard(is.factor) %>% 
  str()
## 'data.frame':    150 obs. of  4 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...

some()とevery()は述語が要素のいずれか、またはすべてでTRUEとなるかどうかを判定する。

x <- list(1:5, letters, list(10))
x %>% 
  some(is_character)
## [1] TRUE
x %>% 
  every(is_vector)
## [1] TRUE

detect()は述語が真となる最初の要素を返す、detect_index()はその位置を返す。

x <- sample(10)
x
##  [1]  9  3  7  8  5  4 10  6  2  1
x %>% 
  detect(~ . > 5)
## [1] 9
x %>% 
  detect_index(~ . > 5)
## [1] 1

head_while()とtail_while()は述語がTRUEである先頭または末尾からの要素列を返す。

x %>% 
  head_while(~ . > 5)
## [1] 9
x %>% 
  tail_while(~ . > 5)
## integer(0)

1.5.3.8 reduceとaccumulate

reduce()は「2項」関数（2つの主入力がある）をとり、単一要素が残るまでリストに繰り返しその関数を適用する。例えばデータフレームのリストに対して、要素を結合して単一データフレームにする場合に使用。

dfs <- list(
  age = tibble(name = "John", age = 30),
  sex = tibble(name = c("John", "Mary"), sex = c("M", "F")),
  trt = tibble(name = "Mary", treatment = "A")
)
dfs %>% reduce(full_join)
## Joining with `by = join_by(name)`
## Joining with `by = join_by(name)`
## # A tibble: 2 × 4
##   name    age sex   treatment
##   <chr> <dbl> <chr> <chr>    
## 1 John     30 M     <NA>     
## 2 Mary     NA F     A

あるいはベクトルのリストがあり、共通部分を求めたい場合。

vs <- list(
  c(1, 3, 5, 6, 10),
  c(1, 2, 3, 7, 8, 10),
  c(1, 2, 3, 4, 8, 9, 10)
)
vs %>% reduce(intersect)  ## intersect()は両方に含まれている要素を抽出する
## [1]  1  3 10

accumulate()もreduce()と同様に機能するが、中間結果を全て保存する。累積和を求めることができる。

x <- sample(10)
x
##  [1]  9  8  1  6  4  2 10  3  7  5
x %>% accumulate(`+`)  ## accumulateは中間結果を全て保存する
##  [1]  9 17 18 24 28 30 40 43 50 55

1.6 文字列の処理

Rにおける文字列の操作と正規表現（regular expression）。stringrパッケージを使用する。

library(tidyverse)
library(stringr)
data("words")

特殊文字はバックスラッシュ\を使ってエスケープする。文字列の中身そのものを表示させるにはwriteLines()。

double_quote <- "\""  ## バックスラッシュでエスケープ
double_quote
## [1] "\""
writeLines(double_quote)  ## 文字列の中身そのものを調べるときはwriteLines
## "
return <- "a \nb "  ## 改行
tab <- "a \t b"  ## タブ
mu <- "\u00b5"
writeLines(return)
## a 
## b
writeLines(tab)
## a     b
writeLines(mu)
## µ
?"'"  ## 特殊文字一覧

1.6.1 stringr関数

文字列を操作するstringrの関数は全て「str_」で始まるので、RStudioの自動補完機能で一覧が表示される。

str_length(c("a", "R for data science", NA))  ## 文字列の長さ
## [1]  1 18 NA
str_c("x", "y") ## 文字列の連結
## [1] "xy"
str_c("x", "y", sep = ", " ) ## sep引数を使って区切り文字を制御
## [1] "x, y"
x <- c("abc", NA)
str_c("|-", x, "-|")  ## NAはそのまま
## [1] "|-abc-|" NA
## "NA"と出力するにはstr_replace_naを使用
str_c("|-", str_replace_na(x), "-|") 
## [1] "|-abc-|" "|-NA-|"

str_c()はベクトル化関数で、要素の少ないベクトルを最長ベクトルの要素数に自動的にリサイクルして合わせる。

str_c("prefix-", c("a", "b", "c"), "-suffix") 
## [1] "prefix-a-suffix" "prefix-b-suffix" "prefix-c-suffix"

引数collapse = ""とすると、文字列のベクトルをまとめて一つの文字列にする。

str_c(c("x", "y", "z"), "a", collapse = "")
## [1] "xayaza"
str_c(c("x", "y", "z"), "a")
## [1] "xa" "ya" "za"
str_c(c("x", "y", "z"), collapse = "")
## [1] "xyz"

1.6.1.1 文字列の一部抽出

str_sub()は文字列が短すぎてもエラーにならず、できるだけ部分文字列を取り出す。

x <- c("Apple", "Banana", "Pear")
str_sub(x, 1, 3)  ## 1〜3文字目を取り出す
## [1] "App" "Ban" "Pea"
str_sub(x, -3, -1)  ## 負数は末尾からの位置
## [1] "ple" "ana" "ear"
# str_sub()の代入形式を使って、文字列を変更することもできる
str_sub(x ,1, 1) <- str_to_lower(str_sub(x, 1, 1))
x
## [1] "apple"  "banana" "pear"

後述の正規表現とstr_extract()を使用することもできる。

x <- c("Apple", "Banana", "Pear")
str_extract(x, "^.{2}") # 先頭の２文字を取り出す
## [1] "Ap" "Ba" "Pe"
str_extract(x, ".{3}$") # 末尾の３文字を取り出す
## [1] "ple" "ana" "ear"

1.6.2 正規表現

1.6.2.1 基本マッチ

str_view()は文字列と正規表現を取り、マッチがどうなるかを示す。最も単純なパターンマッチは厳密なマッチ。改行以外のどんな文字ともマッチする「.」。

x <- c("apple", "banana", "pear")
str_view(x, "an") ## 厳密なマッチ
## [2] │ b<an><an>a
str_view(x, ".a.")  ## .は任意の一文字
## [2] │ <ban>ana
## [3] │ p<ear>

正規表現の「.」ではなく、文字ピリオド「.」をマッチさせるには、エスケープして「\.」とする。ただしR上では文字列のエスケープと認識されてしまい、「\.」のような文字列エスケープは無いためエラーとなる。この問題を回避するため正規表現の\.を作るには、R上ではバックスラッシュを2つ重ねて「\\.」とする。参考：https://opur.club/textbook/2018-5-2/, https://www.jaysong.net/RBook/string.html#%E6%AD%A3%E8%A6%8F%E8%A1%A8%E7%8F%BE

# To create the regular expression, we need \\
dot <- "\\."
# But the expression itself only contains one:
writeLines(dot)
## \.
#> \.
# And this tells R to look for an explicit .
str_view(c("abc", "a.c", "bef"), "a\\.c")
## [2] │ <a.c>

文字通りのバックスラッシュとマッチさせるには、「\\\\」と、4つ重ねる必要がある。

x <- "a\\b"
writeLines(x)
## a\b
#> a\b
str_view(x, "\\\\")
## [1] │ a<\>b
マッチ 正規表現 Rでの表現
任意の数字 \d \\d
空白文字（空白, タブ, 改行） \s \\s
記号を除く文字 \w \\w
単語の境界 \b \\b
ピリオド「.」 \. \\.
バックスラッシュ「\」 \\ \\\\

1.6.2.2 アンカー

  • ^: 文字列の先頭とマッチする
  • $: 文字列の末尾とマッチする
x <- c("apple", "banana", "pear")
str_view(x, "^a")  
## [1] │ <a>pple
str_view(x, "a$") 
## [2] │ banan<a>
x <- c("apple pie", "apple", "apple cake")
str_view(x, "apple")
## [1] │ <apple> pie
## [2] │ <apple>
## [3] │ <apple> cake
## 文字列全体とだけマッチするには^と$の両方でアンカーする
str_view(x, "^apple$")  
## [2] │ <apple>

1.6.2.3 文字のクラスと候補

マッチ 正規表現
a, b, cのどれか1文字 [abc]
a,b,c以外ならどの文字とも [^abc]
この並びで現れる (abc) 
# Look for a literal character that normally has special meaning in a regex
str_view(c("abc", "a.c", "a*c", "a c"), "a[.]c")
## [2] │ <a.c>
str_view(c("abc", "a.c", "a*c", "a c"), ".[*]c")
## [3] │ <a*c>
str_view(c("abc", "a.c", "a*c", "a c"), "a[ ]")
## [4] │ <a >c
# 複数の候補パターンから選ぶには候補指定|を使用する
str_view(c("grey", "gray"), "gr(e|a)y") 
## [1] │ <grey>
## [2] │ <gray>

1.6.2.4 繰り返し

意味 正規表現
0か1 ?
1以上 +
0以上 *
n個 {n}
n個以上 {n,}
n個以下 {,n}
m個以上、n個以下 {m,n} 
x <- "1888 is the longest year in Roman numerals: MDCCCLXXXVIII"
str_view(x, "CC?")
## [1] │ 1888 is the longest year in Roman numerals: MD<CC><C>LXXXVIII
str_view(x, "CC+")
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, 'C[LX]+')　　## [LX]はLXのどれかとマッチする
## [1] │ 1888 is the longest year in Roman numerals: MDCC<CLXXX>VIII
str_view(x, "C{2}")  ## 2回繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CC>CLXXXVIII
str_view(x, "C{2,}")  ## 2回以上繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, "C{2,3}")  ## 2〜3回繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, 'C{2,3}?')  ## defaultでは最長文字列と最長マッチする、?を後ろにつけて最短文字列とマッチさせることもできる
## [1] │ 1888 is the longest year in Roman numerals: MD<CC>CLXXXVIII
str_view(x, 'C[LX]+?')
## [1] │ 1888 is the longest year in Roman numerals: MDCC<CL>XXXVIII

1.6.2.5 正規表現の例

意味 正規表現
先頭がA〜Fで始まる ^[A-F]
末尾が数字の1〜9で終わる [1-9]$
何らかの文字が0回以上出現 .*
3桁の数字 [0-9]{3}
3桁-3桁-4桁の電話番号 ^[0-9]{3}-[0-9]{3}-[0-9]{4}$

1.6.2.6 グループ化と後方参照

()でグループ化すると、キャプチャも同時に行われる。キャプチャとは、()内の文字列を記憶し、後で参照できるようにする機能。記憶した文字列は、後で「\1」, 「\2」, 「\3」,…のように指定することで利用することができる。正規表現内の()の数だけグループができる。以下のwordsは980個の英単語のリスト。

## ()は番号付きのcapturing groupを作成、\1, \2によりそれらを後方参照できる
str_view(fruit, "(..)\\1", match = TRUE)
##  [4] │ b<anan>a
## [20] │ <coco>nut
## [22] │ <cucu>mber
## [41] │ <juju>be
## [56] │ <papa>ya
## [73] │ s<alal> berry
str_view(fruit, "(.)(.)\\2\\1", match = TRUE)
##  [5] │ bell p<eppe>r
## [17] │ chili p<eppe>r
str_view(words, "(.).\\1.\\1", match = TRUE)
## [265] │ <eleve>n

1.6.2.7 マッチの可否

文字ベクトルがパターンにマッチするかどうかの確認は、str_detect()を使う。

x <- c("apple", "banana", "pear")
str_detect(x, "e")  ## 入力と同じ長さの論理ベクトルを返す
## [1]  TRUE FALSE  TRUE

数値表現ではTRUEは1, FALSEは0としてカウントされる。sum()とmean()でそれぞれ、個数と割合がわかる。

str_detect()とgrepl()はほぼ同じ、文法が異なる。

x <- c("apple", "banana", "pear")
grepl("e", x)
## [1]  TRUE FALSE  TRUE
# 一般単語でtで始まる単語の個数
sum(str_detect(words, "^t"))
## [1] 65
# 一般単語で母音で終わる単語の割合
# [aeiou]はa, e, i, o, uのいずれか、$は末尾を指定するアンカー
mean(str_detect(words, "[aeiou]$"))  
## [1] 0.2765306
sum(str_detect(words, "[aeiou]$"))
## [1] 271

母音を含まない単語を探し出す２つの方法。前者の方がわかりやすい。

# 少なくとも母音を1つ含む単語をすべて探し出して、補集合を取る
no_vowels_1 <- !str_detect(words, "[aeiou]")
# 子音だけからなる単語を全て探し出す
# [^aeiou]はaeiou以外、^と$で囲むと文字列全体、+は一回以上の繰り返し
no_vowels_2 <- str_detect(words, "^[^aeiou]+$")
identical(no_vowels_1, no_vowels_2)  ## identicalで比較
## [1] TRUE

1.6.2.8 パターンにマッチする要素の選択

str_detect()を使用してから、論理ベクトルでリストの部分集合を取る。またはstr_subse()を使う。

# 論理ベクトルで部分集合を取る
words[str_detect(words, "x$")]
## [1] "box" "sex" "six" "tax"
# str_subsetを使用
str_subset(words, "x$")  
## [1] "box" "sex" "six" "tax"

通常は対象はデータフレームの列にあるので、filter()を使う。

df <- tibble(
  word = words,
# seq_along(x)はベクトルxの長さと同じ数の数列を作成、1:length(x)と等価、インデックスを簡単に作れる
  i = seq_along(word) 
)
df %>% 
  filter(str_detect(word, "x$"))
## # A tibble: 4 × 2
##   word      i
##   <chr> <int>
## 1 box     108
## 2 sex     747
## 3 six     772
## 4 tax     841

str_count()は単なるマッチの可否ではなく、文字列にマッチがいくつあるかを示す。

x <- c("apple", "banana", "pear")
str_count(x, "a")  ## 文字列にマッチがいくつあるか
## [1] 1 3 1
# 平均すると、1つの単語にどれだけ母音があるか?
mean(str_count(words, "[aeiou]"))
## [1] 1.991837
df %>% 
  mutate(
    vowels = str_count(word, "[aeiou]"),
    consonants = str_count(word, "[^aeiou]")
  )
## # A tibble: 980 × 4
##    word         i vowels consonants
##    <chr>    <int>  <int>      <int>
##  1 a            1      1          0
##  2 able         2      2          2
##  3 about        3      3          2
##  4 absolute     4      4          4
##  5 accept       5      2          4
##  6 account      6      3          4
##  7 achieve      7      4          3
##  8 across       8      2          4
##  9 act          9      1          2
## 10 active      10      3          3
## # ℹ 970 more rows

マッチが重なることは絶対に無いことに注意。

str_count("abababa", "aba")
## [1] 2
str_view_all("abababa", "aba")
## Warning: `str_view_all()` was deprecated in stringr 1.5.0.
## ℹ Please use `str_view()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] │ <aba>b<aba>

1.6.2.9 マッチの抽出

実際にマッチしたテキストを抽出するにはstr_extract()を使う。解析サンプルとしてstringr::sentencesに格納されているHarvard英語コーパスを使用。720個の文がリストに格納されている。

length(sentences)
## [1] 720
head(sentences)
## [1] "The birch canoe slid on the smooth planks." 
## [2] "Glue the sheet to the dark blue background."
## [3] "It's easy to tell the depth of a well."     
## [4] "These days a chicken leg is a rare dish."   
## [5] "Rice is often served in round bowls."       
## [6] "The juice of lemons makes fine punch."

色名を含む文を全て探す。最初に色名ベクトルを作成し、それを一つの正規表現にする。

colors <- c("red", "orange", "yellow", "green", "blue", "purple")
color_match <- str_c(colors, collapse = "|")
color_match
## [1] "red|orange|yellow|green|blue|purple"

色を含む文の部分集合リストを作り、どの色か調べるために色を抽出する。

has_color <- str_subset(sentences, color_match) ## 色を含む文を選択してリストを作る
matches <- str_extract(has_color, color_match)  ## 色を抽出
# head(matches)
str_view_all(has_color, color_match)
##  [1] │ Glue the sheet to the dark <blue> background.
##  [2] │ Two <blue> fish swam in the tank.
##  [3] │ The colt rea<red> and threw the tall rider.
##  [4] │ The wide road shimme<red> in the hot sun.
##  [5] │ See the cat glaring at the sca<red> mouse.
##  [6] │ A wisp of cloud hung in the <blue> air.
##  [7] │ Leaves turn brown and <yellow> in the fall.
##  [8] │ He orde<red> peach pie with ice cream.
##  [9] │ Pure b<red> poodles have curls.
## [10] │ The spot on the blotter was made by <green> ink.
## [11] │ Mud was spatte<red> on the front of his white shirt.
## [12] │ The sofa cushion is <red> and of light weight.
## [13] │ The sky that morning was clear and bright <blue>.
## [14] │ Torn scraps litte<red> the stone floor.
## [15] │ The doctor cu<red> him with these pills.
## [16] │ The new girl was fi<red> today at noon.
## [17] │ The third act was dull and ti<red> the players.
## [18] │ A <blue> crane is a tall wading bird.
## [19] │ Live wires should be kept cove<red>.
## [20] │ It is hard to erase <blue> or <red> ink.
## ... and 37 more

str_extract()は最初のマッチしか抽出しないことに注意。

# マッチを2つ以上含む文を選択
more <- sentences[str_count(sentences, color_match) >1] 
str_view_all(more, color_match)
## [1] │ It is hard to erase <blue> or <red> ink.
## [2] │ The <green> light in the brown box flicke<red>.
## [3] │ The sky in the west is tinged with <orange> <red>.
str_extract(more, color_match)
## [1] "blue"   "green"  "orange"

マッチを全て取得するにはstr_extract_all()を使用する。返り値が2つ以上になるので、ベクトルではなくリストが返ってくる。

str_extract_all(more, color_match)
## [[1]]
## [1] "blue" "red" 
## 
## [[2]]
## [1] "green" "red"  
## 
## [[3]]
## [1] "orange" "red"
# simplify = Tを使えば, str_extract_all()が少ないマッチでも最多マッチに合わせた行列を返す
str_extract_all(more, color_match, simplify = TRUE)
##      [,1]     [,2] 
## [1,] "blue"   "red"
## [2,] "green"  "red"
## [3,] "orange" "red"
x <- c("a", "a b", "a b c")
str_extract_all(x, "[a-z]", simplify = TRUE)
##      [,1] [,2] [,3]
## [1,] "a"  ""   ""  
## [2,] "a"  "b"  ""  
## [3,] "a"  "b"  "c"

色名を単語の途中ではなく、独立した単語として抽出するには、色名を単語の境界を示す\bで挟む。

colors2 <- c("\\bred", "orange", "yellow", "green", "blue", "purple\\b")
color_match2 <- str_c(colors2, collapse = "\\b|\\b")
color_match2
## [1] "\\bred\\b|\\borange\\b|\\byellow\\b|\\bgreen\\b|\\bblue\\b|\\bpurple\\b"
has_color2 <- str_subset(sentences, color_match2)
str_view_all(has_color2, color_match2)
##  [1] │ Glue the sheet to the dark <blue> background.
##  [2] │ Two <blue> fish swam in the tank.
##  [3] │ A wisp of cloud hung in the <blue> air.
##  [4] │ Leaves turn brown and <yellow> in the fall.
##  [5] │ The spot on the blotter was made by <green> ink.
##  [6] │ The sofa cushion is <red> and of light weight.
##  [7] │ The sky that morning was clear and bright <blue>.
##  [8] │ A <blue> crane is a tall wading bird.
##  [9] │ It is hard to erase <blue> or <red> ink.
## [10] │ The lamp shone with a steady <green> flame.
## [11] │ The box is held by a bright <red> snapper.
## [12] │ The houses are built of <red> clay bricks.
## [13] │ The <red> tape bound the smuggled food.
## [14] │ Hedge apples may stain your hands <green>.
## [15] │ The plant grew large and <green> in the window.
## [16] │ The <purple> tie was ten years old.
## [17] │ Bathe and relax in the cool <green> grass.
## [18] │ The lake sparkled in the <red> hot sun.
## [19] │ Mark the spot with a sign painted <red>.
## [20] │ The couch cover and hall drapes were <blue>.
## ... and 9 more

「ing」で終わる単語を抽出する。\\bで単語の区切りを設定。\w*は任意の文字を任意の回数繰り返す。その後にingが続く。

end_ing = "\\b\\w*ing\\b"
# end_ing = "\\b[^ ]+ing\\b" # 別法："[^ ]+"は「空白以外の文字の1つ以上の連なり」
ing_match <- str_subset(sentences, end_ing)
# str_view_all(ing_match, end_ing)
head(str_extract_all(ing_match, end_ing))
## [[1]]
## [1] "spring"
## 
## [[2]]
## [1] "evening"
## 
## [[3]]
## [1] "morning"
## 
## [[4]]
## [1] "winding"
## 
## [[5]]
## [1] "living"
## 
## [[6]]
## [1] "king"

1.6.2.10 グループ化したマッチ

文章から名詞を抽出するために、“a”または”the”の後に来る全ての単語を探す。()を使うことにより、2つのマッチがグループになる。

## a または the の後に、空白以外の文字の一つ以上の連なりを名詞の近似とする
noun <- "(a|the) ([^ ]+)"  
has_noun <- sentences %>% str_subset(noun) %>% head(10)
has_noun %>% str_extract(noun)
##  [1] "the smooth" "the sheet"  "the depth"  "a chicken"  "the parked"
##  [6] "the sun"    "the huge"   "the ball"   "the woman"  "a helps"

str_extract()は完全マッチ、str_match()は個別の部分マッチを示す、行列を返し第1列が完全マッチ、第2列以降がそれぞれのグループstr_extract()と同様、各文字列全てにマッチしたいなら、str_match_all()が必要。

has_noun %>% str_match(noun)
##       [,1]         [,2]  [,3]     
##  [1,] "the smooth" "the" "smooth" 
##  [2,] "the sheet"  "the" "sheet"  
##  [3,] "the depth"  "the" "depth"  
##  [4,] "a chicken"  "a"   "chicken"
##  [5,] "the parked" "the" "parked" 
##  [6,] "the sun"    "the" "sun"    
##  [7,] "the huge"   "the" "huge"   
##  [8,] "the ball"   "the" "ball"   
##  [9,] "the woman"  "the" "woman"  
## [10,] "a helps"    "a"   "helps"

データフレームを検索する場合はtidyr::extractを使用、str_extract()と同様に作用するが、マッチに名前を付ける必要がある。

tibble(sentence = sentences) %>% 
  tidyr::extract(
    sentence, c("article", "noun"), "(a|the) ([^ ]+)", 
    remove = FALSE
  )
## # A tibble: 720 × 3
##    sentence                                    article noun   
##    <chr>                                       <chr>   <chr>  
##  1 The birch canoe slid on the smooth planks.  the     smooth 
##  2 Glue the sheet to the dark blue background. the     sheet  
##  3 It's easy to tell the depth of a well.      the     depth  
##  4 These days a chicken leg is a rare dish.    a       chicken
##  5 Rice is often served in round bowls.        <NA>    <NA>   
##  6 The juice of lemons makes fine punch.       <NA>    <NA>   
##  7 The box was thrown beside the parked truck. the     parked 
##  8 The hogs were fed chopped corn and garbage. <NA>    <NA>   
##  9 Four hours of steady work faced us.         <NA>    <NA>   
## 10 A large size in stockings is hard to sell.  <NA>    <NA>   
## # ℹ 710 more rows

1.6.2.11 マッチの置換

x <- c("apple", "pear", "banana")
str_replace(x, "[aeiou]", "-")　　## 最初のマッチを置換
## [1] "-pple"  "p-ar"   "b-nana"
str_replace_all(x, "[aeiou]", "-")　## 全てのマッチを置換
## [1] "-ppl-"  "p--r"   "b-n-n-"
x <- c("1 house", "2 cars", "3 people")
str_replace_all(x, c("1" = "one", "2" = "two", "3" = "three"))  ## 名前付きベクトルを指定すると複数の置換ができる
## [1] "one house"    "two cars"     "three people"

固定文字列で置換する代わりに、後方参照を使ってマッチした要素を挿入することもできる。

sentences %>% 
  str_replace("([^ ]+) ([^ ]+) ([^ ]+)", "\\1 \\3 \\2") %>% ## 2番目と3番目の単語の順序を入れ替える
  head(5)
## [1] "The canoe birch slid on the smooth planks." 
## [2] "Glue sheet the to the dark blue background."
## [3] "It's to easy tell the depth of a well."     
## [4] "These a days chicken leg is a rare dish."   
## [5] "Rice often is served in round bowls."

1.6.2.12 分割

str_split()を使うと文字列を分割できる。空白で分割すると、文を単語に分割できる。文ごとに要素の数が異なるのでリストを返す。

sentences %>% head(5) %>% str_split(" ")
## [[1]]
## [1] "The"     "birch"   "canoe"   "slid"    "on"      "the"     "smooth" 
## [8] "planks."
## 
## [[2]]
## [1] "Glue"        "the"         "sheet"       "to"          "the"        
## [6] "dark"        "blue"        "background."
## 
## [[3]]
## [1] "It's"  "easy"  "to"    "tell"  "the"   "depth" "of"    "a"     "well."
## 
## [[4]]
## [1] "These"   "days"    "a"       "chicken" "leg"     "is"      "a"      
## [8] "rare"    "dish."  
## 
## [[5]]
## [1] "Rice"   "is"     "often"  "served" "in"     "round"  "bowls."
"a|b|c|d" %>% 
  str_split("\\|") %>% ## 長さ1のベクトルで返したいなら、リストの最初の要素[[1]]だけを抽出
  .[[1]]
## [1] "a" "b" "c" "d"
sentences %>% ## simplity =TRUEで行列を返す
  head(5) %>% 
  str_split(" ", simplify = TRUE)
##      [,1]    [,2]    [,3]    [,4]      [,5]  [,6]    [,7]     [,8]         
## [1,] "The"   "birch" "canoe" "slid"    "on"  "the"   "smooth" "planks."    
## [2,] "Glue"  "the"   "sheet" "to"      "the" "dark"  "blue"   "background."
## [3,] "It's"  "easy"  "to"    "tell"    "the" "depth" "of"     "a"          
## [4,] "These" "days"  "a"     "chicken" "leg" "is"    "a"      "rare"       
## [5,] "Rice"  "is"    "often" "served"  "in"  "round" "bowls." ""           
##      [,9]   
## [1,] ""     
## [2,] ""     
## [3,] "well."
## [4,] "dish."
## [5,] ""

要素の最大個数を指定できる

fields <- c("Name: Hadley", "Country: NZ", "Age: 35")
fields %>% str_split(": ", n = 2, simplify = TRUE)
##      [,1]      [,2]    
## [1,] "Name"    "Hadley"
## [2,] "Country" "NZ"    
## [3,] "Age"     "35"
x <- "This is a sentence.  This is another sentence."
str_view_all(x, boundary("word")) ## パターンではなく、文字、行、文、後のboundary()で分割できる
## [1] │ <This> <is> <a> <sentence>.  <This> <is> <another> <sentence>.
str_split(x, " ")[[1]]
## [1] "This"      "is"        "a"         "sentence." ""          "This"     
## [7] "is"        "another"   "sentence."
str_split(x, boundary("word"))[[1]]
## [1] "This"     "is"       "a"        "sentence" "This"     "is"       "another" 
## [8] "sentence"

1.6.2.13 他の種類のパターン

文字列のパターンを使えば自動的にregex()にラップされる。

str_view(fruit, "nana")
## [4] │ ba<nana>
# Is shorthand for
str_view(fruit, regex("nana"))
## [4] │ ba<nana>

regex()の他の引数を使ってマッチの詳細を制御できる。

bananas <- c("banana", "Banana", "BANANA")
str_view(bananas, "banana")
## [1] │ <banana>
str_view(bananas, regex("banana", ignore_case = TRUE))  ## ignore_caseで大文字小文字を無視
## [1] │ <banana>
## [2] │ <Banana>
## [3] │ <BANANA>

multilineは^と$が文字列全体の先頭と末尾ではなく、各行の先頭と末尾にマッチする。

x <- "Line 1\nLine 2\nLine 3"
str_extract_all(x, "^Line")[[1]]
## [1] "Line"
str_extract_all(x, regex("^Line", multiline = TRUE))[[1]] 
## [1] "Line" "Line" "Line"

commentsは正規表現を読みやすくするためにコメントと空白を使う。

phone <- regex("
  \\(?     # optional opening parens
  (\\d{3}) # area code
  [) -]?   # optional closing parens, space, or dash
  (\\d{3}) # another three numbers
  [ -]?    # optional space or dash
  (\\d{3}) # three more numbers
  ", comments = TRUE) 
str_match("514-791-8141", phone)
##      [,1]          [,2]  [,3]  [,4] 
## [1,] "514-791-814" "514" "791" "814"

2 データの読み込みと保存

データ解析の対象となるような大きなデータセットは、通常は手作業で入力することは無く、ファイルを読み込んでデータフレームオブジェクトとして保存する。

2.1 RStudioによるデータの読み込み

Environmental paneのImport DatasetボタンからcsvファイルやExcelファイルを読み込むことができる。

2.2 関数によるデータの読み込みと保存

read.csv()やread_excel()関数を使用する。

# library(readxl) ## excelやcsv fileを読み込む
# library(openxlsx)   ## excelファイルを読み書きする
seiseki <- read.csv("seiseki.csv", header =T, sep = ",", stringsAsFactors = F)
test_data <- read_excel("test_data.xlsx", sheet = "200302")  # 読み込むシートも指定できる
head(seiseki)
##   kokugo shakai sugaku rika ongaku bijutu taiiku gika eigo
## 1     30     43     51   63     60     66     37   44   20
## 2     39     21     49   56     70     72     56   63   16
## 3     29     30     23   57     69     76     33   54    6
## 4     95     87     77  100     77     82     78   96   87
## 5     70     71     78   67     72     82     46   63   44
## 6     67     53     56   61     61     76     70   66   40
head(test_data)
## # A tibble: 6 × 8
##   Date       Filename      `# of z-slice` defocus `cells with granule`
##   <chr>      <chr>                  <dbl>   <dbl>                <dbl>
## 1 2020-03-02 ade12_006.nd2             22      NA                  213
## 2 2020-03-02 ade12_007.nd2             22      NA                  258
## 3 2020-03-02 ade12_008.nd2             22      NA                  254
## 4 2020-03-02 ade12_020.nd2             22       1                   NA
## 5 2020-03-02 ade12_021.nd2             22      NA                  307
## 6 2020-03-02 ade12_022.nd2             22      NA                  308
## # ℹ 3 more variables: `cell without granule` <dbl>, `Total cell` <dbl>,
## #   pct_cell_with_granules <dbl>

Excelファイルやテキストファイル上で該当するデータ範囲をクリップボードにコピーして、データフレームとして読み込みことができる。 WindowsとMacでコマンドが異なる。

# clip_tab <- read.table(pipe("pbpaste"), sep = "", header = T)  # Macの場合, sepは区切り文字を指定、header = Tは先頭行を列名にする
# clip_tab <- read.table("clipboard", sep = "", header = T)  # Windowsの場合

Rのデータフレームをcsv形式で保存するにはwrite.csv()を、Excel形式で保存するにはwrite.xlsx()をそれぞれ使用。

write.csv(test_data, file = "test_data_csv.csv")
write.xlsx(seiseki, file = "seiseki.xlsx", colNames = T, sheetName = "seiseki")

2.3 複数のデータフレームの結合

2つ以上のデータフレームを縦に連結するにはbind_rows()を使用する。同じ名前の列を結合する。

df1 <- read_excel("test_data.xlsx", sheet = "200302")
str(df1)
## tibble [30 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Date                  : chr [1:30] "2020-03-02" "2020-03-02" "2020-03-02" "2020-03-02" ...
##  $ Filename              : chr [1:30] "ade12_006.nd2" "ade12_007.nd2" "ade12_008.nd2" "ade12_020.nd2" ...
##  $ # of z-slice          : num [1:30] 22 22 22 22 22 22 22 22 22 22 ...
##  $ defocus               : num [1:30] NA NA NA 1 NA NA NA NA NA NA ...
##  $ cells with granule    : num [1:30] 213 258 254 NA 307 308 260 94 144 110 ...
##  $ cell without granule  : num [1:30] 75 91 109 NA 70 80 55 144 224 181 ...
##  $ Total cell            : num [1:30] 288 349 363 NA 377 388 315 238 368 291 ...
##  $ pct_cell_with_granules: num [1:30] 74 73.9 70 NA 81.4 ...
df2 <- read_excel("test_data.xlsx", sheet = "200304")
str(df2)
## tibble [18 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Date                  : chr [1:18] "2020-03-04" "2020-03-04" "2020-03-04" "2020-03-04" ...
##  $ Filename              : chr [1:18] "ade12_005.nd2" "ade12_006.nd2" "ade12_007.nd2" "ade12_hyp_008.nd2" ...
##  $ # of z-slice          : num [1:18] 22 22 22 22 22 22 22 22 22 22 ...
##  $ defocus               : num [1:18] NA NA NA NA NA NA NA NA NA NA ...
##  $ cells with granule    : num [1:18] 235 249 192 154 182 146 194 170 120 179 ...
##  $ cell without granule  : num [1:18] 53 56 57 56 52 52 66 45 50 41 ...
##  $ Total cell            : num [1:18] 288 305 249 210 234 198 260 215 170 220 ...
##  $ pct_cell_with_granules: num [1:18] 81.6 81.6 77.1 73.3 77.8 ...
df_merge <- bind_rows(df1, df2)
str(df_merge)
## tibble [48 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Date                  : chr [1:48] "2020-03-02" "2020-03-02" "2020-03-02" "2020-03-02" ...
##  $ Filename              : chr [1:48] "ade12_006.nd2" "ade12_007.nd2" "ade12_008.nd2" "ade12_020.nd2" ...
##  $ # of z-slice          : num [1:48] 22 22 22 22 22 22 22 22 22 22 ...
##  $ defocus               : num [1:48] NA NA NA 1 NA NA NA NA NA NA ...
##  $ cells with granule    : num [1:48] 213 258 254 NA 307 308 260 94 144 110 ...
##  $ cell without granule  : num [1:48] 75 91 109 NA 70 80 55 144 224 181 ...
##  $ Total cell            : num [1:48] 288 349 363 NA 377 388 315 238 368 291 ...
##  $ pct_cell_with_granules: num [1:48] 74 73.9 70 NA 81.4 ...
df_list <- list(df1, df2)
df_merge2 <- bind_rows(df_list) # bind_rowsにわたす引数はリストでも可

2.4 多数のデータファイルの読み込み

2.4.1 単一のデータフレームにまとめる

多数のデータファイルからデータを読み込んで、1つのデータフレームに集約させる場合は次の3段階で行う：

  1. 読み込むcsv/Excelファイルをワーキングディレクトリの下のデータフォルダにまとめる、サブフォルダに分かれていても良い。
  2. データフォルダ中のデータファイルのパスをサブフォルダも含めて取得してリストを作る。
  3. リストを使用してデータファイルを一括で読み込み、データフレーム化した後に縦に結合する。

以下の例ではデータフォルダmaxima_dataの下に5つのサブフォルダがあり、それぞれに複数のcsvファイルが格納されている。
Excelファイルを読み込む場合は.csvを.xlsxに変える。

data_dir <- "maxima_data"  # データフォルダの指定
# list.files()関数でフォルダ内のファイル名リストを取得、recursiveでサブフォルダも探索する、拡張子も指定できる
csv_list <- list.files(data_dir, full.names = T, recursive = T, pattern = "*.csv") 
str(csv_list) # 合計で84個のcsvファイル
##  chr [1:84] "maxima_data/230222_maxima_counts/maxima_counts_Ade4-mNG_sc-lf_003_prom100.csv" ...
apply_res <- csv_list %>%
  lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) # lapplyはリストの各要素に対して演算を行い、結果をリストで返す
length(apply_res) # 84個のデータフレームが格納されている
## [1] 84
csv_merge <- bind_rows(apply_res)  # bind_rowsでapply_res中のデータフレームを連結
str(csv_merge)
## 'data.frame':    1260 obs. of  10 variables:
##  $ X             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Count         : int  29 26 29 39 45 58 69 80 83 101 ...
##  $ total_cell_num: int  90 90 90 90 90 90 90 90 90 90 ...
##  $ count_per_cell: num  0.322 0.289 0.322 0.433 0.5 ...
##  $ prominence    : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ BC_Mean       : num  17.3 17.3 17.3 17.3 17.3 ...
##  $ BC_Mode       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ BC_Median     : int  14 14 14 14 14 14 14 14 14 14 ...
##  $ date          : int  230222 230222 230222 230222 230222 230222 230222 230222 230222 230222 ...
##  $ file          : chr  "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" ...
# apply_resを作成せずに一気に連結する場合
csv_merge2 <- csv_list %>%
  lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) %>% 
  bind_rows()
str(csv_merge2)
## 'data.frame':    1260 obs. of  10 variables:
##  $ X             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Count         : int  29 26 29 39 45 58 69 80 83 101 ...
##  $ total_cell_num: int  90 90 90 90 90 90 90 90 90 90 ...
##  $ count_per_cell: num  0.322 0.289 0.322 0.433 0.5 ...
##  $ prominence    : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ BC_Mean       : num  17.3 17.3 17.3 17.3 17.3 ...
##  $ BC_Mode       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ BC_Median     : int  14 14 14 14 14 14 14 14 14 14 ...
##  $ date          : int  230222 230222 230222 230222 230222 230222 230222 230222 230222 230222 ...
##  $ file          : chr  "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" ...

2.4.2 単一のExcelファイルの各シートに出力

読み込んだ多数のデータファイル（.csv or .xlsx）をデータフレームに変換して、1つのExcelファイルのワークシートにそれぞれ出力する場合。

data_dir2 <- "maxima_data/230228_maxima_counts"
csv_list2 <- list.files(data_dir2, full.names = T, recursive = T, pattern = "*.csv") 
sheetname_list <- basename(csv_list2) %>%  # basenameは上位のフォルダ名を含まないファイル名を返す
  str_replace_all(c("maxima_counts_" = "", ".csv" = "")) %>% # str_replace_all()で複数マッチを一度に置換、接頭語と拡張子を除く
  str_sub (., start = 1, end = 31) # シート名は31文字以内の制限があるので、さらに先頭から31文字を切り出す
str(sheetname_list)
##  chr [1:24] "Ade4-mNG_sc-lf_003_prom100" "Ade4-mNG_sc-lf_003_prom110" ...
names(csv_list2) <- sheetname_list  # csv_list2の要素の名前を、短縮したファイル名のリストsheetname_listで置換

file_list <- csv_list2 %>%
  lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) # bind_rowsで結合せずにリストにデータフレームを格納
length(file_list)
## [1] 24
file_list %>% write.xlsx(., file = "assembled.xlsx", colNames = T)  ## xlsxファイルに書き出し、各シートが各ファイルに対応

3 データの整理と変換

3.1 用語の定義

  • 変数 variable：測定できる量、質、または特性。値が離散的なカテゴリ変数と、値が連続的な数値変数がある。
  • 値 value：測定時の変数の状態。 変数の値は測定ごとに変化する場合がある。
  • 観察/観測 observation：同様の条件下で行われた一連の測定。通常、観察ですべての測定を同時に同じ対象に対して行う。 観測には複数の値が含まれ、それぞれ異なる変数に関連付けられている。 観測をデータポイントと呼ぶことがある。 表形式のデータは一連の値であり、それぞれが変数と観察に関連付けられている。
  • 一意識別子 unique identifier：データ内の全ての行を一意に識別するための変数または変数の組み合わせ。 単一の変数で一意性を表す場合は、基本的に機械的につけた番号やコードになる。
  • id変数 id variable：あるカテゴリについて、個別の要素を識別するための固有の数字を表す変数（例：運転免許証番号、クレジットカード番号など）。

3.2 整理データ

整理データ tidy dataは以下の条件を満たすデータ：

  1. 各変数（variable）に専用の列に配置されており、列名は変数名
  2. 各観測（observation）に専用の行に配置されている
  3. 各値（value）は専用のセルに配置されている
    パッケージtidyverseにはtibble形式の整理データを変換するための優れた関数が含まれている。

3.3 パイプ演算子

パイプ演算子 %>% は左辺の計算結果のオブジェクトを右辺の第1引数に渡す。第1引数以外で呼び出す時は .(ドット演算子)で行う。

# library(magrittr)
3 %>% c(., 3 + .) 
## [1] 3 6
3 %>% c(3 + .) # 右辺の第一引数は省略可能
## [1] 3 6
3 %>% { c(3 + .) }  # 中括弧はパイプが関数内の最初の引数を使用しないようにする 
## [1] 6
mtcars %>% {cor(.$mpg, .$cyl)}
## [1] -0.852162
1:10 %>% mean(.) 
## [1] 5.5
1:10 %>% mean # 右辺の第一引数は省略可能
## [1] 5.5
1:10 %>% plot(x = seq(10, 100, 10), y = .)

Dollar演算子 %$% はデータフレームの列にアクセスできる

iris %>% 
  filter(Species != "setosa") %$% 
  lm(Sepal.Length ~ Sepal.Width) %>%   ## lm(.$Sepal.Length ~ .$Sepal.Width)と同等
  summary()
## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0032 -0.3877 -0.0774  0.3200  1.7381 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.0934     0.4844   6.387 5.70e-09 ***
## Sepal.Width   1.1033     0.1675   6.585 2.27e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5547 on 98 degrees of freedom
## Multiple R-squared:  0.3068, Adjusted R-squared:  0.2997 
## F-statistic: 43.36 on 1 and 98 DF,  p-value: 2.27e-09

3.4 データの整理

一般的なデータセットは整理データの形式になっていないことの方が多い。
データ解析の前にデータを整理する必要がある。まず何が変数で、何が観測かを見極めることが重要。

3.4.1 KeyとValueの関係

日付とその日の天気のような、2つの情報の組を考える。個々の組を識別する情報をKey「キー」と呼び、もう1つをValue「値」と呼ぶ。 Key-Value型データベースは、このようなkeyとvalueの組を保存するためのデータベース、辞書データやハッシュテーブルとしても知られる。 複数のkeyでvalueを識別する場合もある。

3.4.2 ワイド形式とロング形式

列方向すなわち横方向に伸びるようにまとめられたデータをwide format（ワイド形式）、 行方向すなわち縦方向に表が伸びていくようにまとめられたデータをlong format（ロング形式）という。ワイドかロングかは相対的な関係である。 ワイド形式は多数のkeyが列名になっており、valueを示す列は存在しないことが多い。ロング形式は少数のkeyと1つのvalueの列から構成されることが多い。

library(tidyverse)
table4a  # tidyverseに収録されているサンプルデータ、ワイド形式、変数の「値」が列名になっている
## # A tibble: 3 × 3
##   country     `1999` `2000`
##   <chr>        <dbl>  <dbl>
## 1 Afghanistan    745   2666
## 2 Brazil       37737  80488
## 3 China       212258 213766
table4a %>% gather("1999":"2000", key =  "year", value = "cases")  # gatherでロング形式に変換、変換する列、列名にする変数名key、列名にする値名valueをそれぞれ指定
## # A tibble: 6 × 3
##   country     year   cases
##   <chr>       <chr>  <dbl>
## 1 Afghanistan 1999     745
## 2 Brazil      1999   37737
## 3 China       1999  212258
## 4 Afghanistan 2000    2666
## 5 Brazil      2000   80488
## 6 China       2000  213766
table4a %>% pivot_longer(cols = c("1999", "2000"), names_to = "year", values_to = "cases") 　#　同機能のより新しい関数 
## # A tibble: 6 × 3
##   country     year   cases
##   <chr>       <chr>  <dbl>
## 1 Afghanistan 1999     745
## 2 Afghanistan 2000    2666
## 3 Brazil      1999   37737
## 4 Brazil      2000   80488
## 5 China       1999  212258
## 6 China       2000  213766
table2  # 列typeに変数casesとpopulationが入っており、変数の値は列countに入っている
## # A tibble: 12 × 4
##    country      year type            count
##    <chr>       <dbl> <chr>           <dbl>
##  1 Afghanistan  1999 cases             745
##  2 Afghanistan  1999 population   19987071
##  3 Afghanistan  2000 cases            2666
##  4 Afghanistan  2000 population   20595360
##  5 Brazil       1999 cases           37737
##  6 Brazil       1999 population  172006362
##  7 Brazil       2000 cases           80488
##  8 Brazil       2000 population  174504898
##  9 China        1999 cases          212258
## 10 China        1999 population 1272915272
## 11 China        2000 cases          213766
## 12 China        2000 population 1280428583
table2 %>% spread(key = type, value = count) # spreadで広げる、慣例として既存の列は""で囲まず、新規に作成する際は""で囲む
## # A tibble: 6 × 4
##   country      year  cases population
##   <chr>       <dbl>  <dbl>      <dbl>
## 1 Afghanistan  1999    745   19987071
## 2 Afghanistan  2000   2666   20595360
## 3 Brazil       1999  37737  172006362
## 4 Brazil       2000  80488  174504898
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583
table2 %>% pivot_wider(names_from = type, values_from = count)  # 同機能のより新しい関数
## # A tibble: 6 × 4
##   country      year  cases population
##   <chr>       <dbl>  <dbl>      <dbl>
## 1 Afghanistan  1999    745   19987071
## 2 Afghanistan  2000   2666   20595360
## 3 Brazil       1999  37737  172006362
## 4 Brazil       2000  80488  174504898
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583

3.4.3 分割と結合

１つの列に２つの変数が含まれている場合、separate()で分ける

table3
## # A tibble: 6 × 3
##   country      year rate             
##   <chr>       <dbl> <chr>            
## 1 Afghanistan  1999 745/19987071     
## 2 Afghanistan  2000 2666/20595360    
## 3 Brazil       1999 37737/172006362  
## 4 Brazil       2000 80488/174504898  
## 5 China        1999 212258/1272915272
## 6 China        2000 213766/1280428583
table3 %>% separate(rate, into = c("cases", "population"))  ## separate()はdefaultで非英数文字で値を分割する
## # A tibble: 6 × 4
##   country      year cases  population
##   <chr>       <dbl> <chr>  <chr>     
## 1 Afghanistan  1999 745    19987071  
## 2 Afghanistan  2000 2666   20595360  
## 3 Brazil       1999 37737  172006362 
## 4 Brazil       2000 80488  174504898 
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583
table3 %>% separate(rate, into = c("cases", "population"), sep = "/") ## sep引数で指定することも可能
## # A tibble: 6 × 4
##   country      year cases  population
##   <chr>       <dbl> <chr>  <chr>     
## 1 Afghanistan  1999 745    19987071  
## 2 Afghanistan  2000 2666   20595360  
## 3 Brazil       1999 37737  172006362 
## 4 Brazil       2000 80488  174504898 
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583
table3 %>% separate(rate, into = c("cases", "population"), convert = T)  ## convertにより数字がより良いint型に変換される
## # A tibble: 6 × 4
##   country      year  cases population
##   <chr>       <dbl>  <int>      <int>
## 1 Afghanistan  1999    745   19987071
## 2 Afghanistan  2000   2666   20595360
## 3 Brazil       1999  37737  172006362
## 4 Brazil       2000  80488  174504898
## 5 China        1999 212258 1272915272
## 6 China        2000 213766 1280428583
table3 %>% separate(year, into = c("century", "year"), sep = 2)  ## sepに整数値を渡すと、文字列の左端を+1, 右端を-1としてその整数の位置で分割する
## # A tibble: 6 × 4
##   country     century year  rate             
##   <chr>       <chr>   <chr> <chr>            
## 1 Afghanistan 19      99    745/19987071     
## 2 Afghanistan 20      00    2666/20595360    
## 3 Brazil      19      99    37737/172006362  
## 4 Brazil      20      00    80488/174504898  
## 5 China       19      99    212258/1272915272
## 6 China       20      00    213766/1280428583

unite()はseparate()の逆の変換。

table5
## # A tibble: 6 × 4
##   country     century year  rate             
##   <chr>       <chr>   <chr> <chr>            
## 1 Afghanistan 19      99    745/19987071     
## 2 Afghanistan 20      00    2666/20595360    
## 3 Brazil      19      99    37737/172006362  
## 4 Brazil      20      00    80488/174504898  
## 5 China       19      99    212258/1272915272
## 6 China       20      00    213766/1280428583
table5 %>% unite("new", century, year, sep = "") ## sepを指定しないデフォルトでは「_」で値が繋がれる
## # A tibble: 6 × 3
##   country     new   rate             
##   <chr>       <chr> <chr>            
## 1 Afghanistan 1999  745/19987071     
## 2 Afghanistan 2000  2666/20595360    
## 3 Brazil      1999  37737/172006362  
## 4 Brazil      2000  80488/174504898  
## 5 China       1999  212258/1272915272
## 6 China       2000  213766/1280428583

3.5 データの変換

dplyrの基本関数：

  • filter()：値から観測値を選び出す
  • arrange()：行を並び替える
  • select()：名前で変数(列)を選ぶ
  • mutate()：既存の変数の関数で新たな変数を作る
  • summarise()：多数の値から単一の要約量を作る

3.5.1 filter

library(tidyverse)
library(nycflights13)
data("flights")
filter(flights, month==1, day==1)  ## 1月1日のフライト、複数の条件は＆になるので両方満たした行が取り出される
## # A tibble: 842 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # ℹ 832 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
filter(flights, month == 11 | month ==12)  ## 11月または12月のフライト
## # A tibble: 55,403 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013    11     1        5           2359         6      352            345
##  2  2013    11     1       35           2250       105      123           2356
##  3  2013    11     1      455            500        -5      641            651
##  4  2013    11     1      539            545        -6      856            827
##  5  2013    11     1      542            545        -3      831            855
##  6  2013    11     1      549            600       -11      912            923
##  7  2013    11     1      550            600       -10      705            659
##  8  2013    11     1      554            600        -6      659            701
##  9  2013    11     1      554            600        -6      826            827
## 10  2013    11     1      554            600        -6      749            751
## # ℹ 55,393 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
# nov_dec <- filter(flights, month %in% c(11, 12))  ## このようにも書ける
filter(flights, arr_delay <= 120, dep_delay <= 120)  ## 出発または到着で120分を超える遅延がなかったフライト
## # A tibble: 316,050 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # ℹ 316,040 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
# filter(flights, !(arr_delay > 120 | dep_delay > 120)) # このようにも書ける
filter(flights, between(dep_time, 0, 6)) ## betweenは範囲を指定、指定した数値も含む
## # A tibble: 155 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     9        2           2359         3      432            444
##  2  2013     1    10        3           2359         4      426            437
##  3  2013     1    13        1           2249        72      108           2357
##  4  2013     1    13        2           2359         3      502            444
##  5  2013     1    13        3           2030       213      340           2350
##  6  2013     1    16        2           2125       157      119           2250
##  7  2013     1    16        3           1946       257      212           2154
##  8  2013     1    22        5           2359         6      452            444
##  9  2013     1    30        3           2159       124      100           2306
## 10  2013     1    31        1           2100       181      124           2225
## # ℹ 145 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>

3.5.2 arrange

# 順序付ける列名の集合を引数に取る。
# 複数の列名を与えると、前列の順序で同じ順序になったものをさらに順序付けるために後列を使用。
arrange(flights, year, month, day)
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
arrange(flights, desc(arr_delay))  ## 降順にするにはdescを使用、NAは常に一番最後になる
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     9      641            900      1301     1242           1530
##  2  2013     6    15     1432           1935      1137     1607           2120
##  3  2013     1    10     1121           1635      1126     1239           1810
##  4  2013     9    20     1139           1845      1014     1457           2210
##  5  2013     7    22      845           1600      1005     1044           1815
##  6  2013     4    10     1100           1900       960     1342           2211
##  7  2013     3    17     2321            810       911      135           1020
##  8  2013     7    22     2257            759       898      121           1026
##  9  2013    12     5      756           1700       896     1058           2020
## 10  2013     5     3     1133           2055       878     1250           2215
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>

3.5.3 select

select(flights, year, month, day) ## 列を名前で選ぶ
## # A tibble: 336,776 × 3
##     year month   day
##    <int> <int> <int>
##  1  2013     1     1
##  2  2013     1     1
##  3  2013     1     1
##  4  2013     1     1
##  5  2013     1     1
##  6  2013     1     1
##  7  2013     1     1
##  8  2013     1     1
##  9  2013     1     1
## 10  2013     1     1
## # ℹ 336,766 more rows
select(flights, year:day)  ## yearとdayの間にある全ての列を選ぶ
## # A tibble: 336,776 × 3
##     year month   day
##    <int> <int> <int>
##  1  2013     1     1
##  2  2013     1     1
##  3  2013     1     1
##  4  2013     1     1
##  5  2013     1     1
##  6  2013     1     1
##  7  2013     1     1
##  8  2013     1     1
##  9  2013     1     1
## 10  2013     1     1
## # ℹ 336,766 more rows
select(flights, -(year:day)) ## yearとdayの間の列以外の全ての列
## # A tibble: 336,776 × 16
##    dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay carrier
##       <int>          <int>     <dbl>    <int>          <int>     <dbl> <chr>  
##  1      517            515         2      830            819        11 UA     
##  2      533            529         4      850            830        20 UA     
##  3      542            540         2      923            850        33 AA     
##  4      544            545        -1     1004           1022       -18 B6     
##  5      554            600        -6      812            837       -25 DL     
##  6      554            558        -4      740            728        12 UA     
##  7      555            600        -5      913            854        19 B6     
##  8      557            600        -3      709            723       -14 EV     
##  9      557            600        -3      838            846        -8 B6     
## 10      558            600        -2      753            745         8 AA     
## # ℹ 336,766 more rows
## # ℹ 9 more variables: flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## #   air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>

select()ではヘルパー関数 start_with()、ends_with()、contains()、num_range()等が使用できる。

select(flights, starts_with("dep"))  # "dep"で始まる列
## # A tibble: 336,776 × 2
##    dep_time dep_delay
##       <int>     <dbl>
##  1      517         2
##  2      533         4
##  3      542         2
##  4      544        -1
##  5      554        -6
##  6      554        -4
##  7      555        -5
##  8      557        -3
##  9      557        -3
## 10      558        -2
## # ℹ 336,766 more rows
select(flights, ends_with("time"))  # "time"で終わる列
## # A tibble: 336,776 × 5
##    dep_time sched_dep_time arr_time sched_arr_time air_time
##       <int>          <int>    <int>          <int>    <dbl>
##  1      517            515      830            819      227
##  2      533            529      850            830      227
##  3      542            540      923            850      160
##  4      544            545     1004           1022      183
##  5      554            600      812            837      116
##  6      554            558      740            728      150
##  7      555            600      913            854      158
##  8      557            600      709            723       53
##  9      557            600      838            846      140
## 10      558            600      753            745      138
## # ℹ 336,766 more rows
select(flights, contains(c("air", "sched", "arr")))  # "air", "sched", "arr"のいずれかを含む
## # A tibble: 336,776 × 6
##    air_time sched_dep_time sched_arr_time arr_time arr_delay carrier
##       <dbl>          <int>          <int>    <int>     <dbl> <chr>  
##  1      227            515            819      830        11 UA     
##  2      227            529            830      850        20 UA     
##  3      160            540            850      923        33 AA     
##  4      183            545           1022     1004       -18 B6     
##  5      116            600            837      812       -25 DL     
##  6      150            558            728      740        12 UA     
##  7      158            600            854      913        19 B6     
##  8       53            600            723      709       -14 EV     
##  9      140            600            846      838        -8 B6     
## 10      138            600            745      753         8 AA     
## # ℹ 336,766 more rows
select(flights, matches("time$"))  # 正規表現にマッチする列名、ここでは"time"で終わる列
## # A tibble: 336,776 × 5
##    dep_time sched_dep_time arr_time sched_arr_time air_time
##       <int>          <int>    <int>          <int>    <dbl>
##  1      517            515      830            819      227
##  2      533            529      850            830      227
##  3      542            540      923            850      160
##  4      544            545     1004           1022      183
##  5      554            600      812            837      116
##  6      554            558      740            728      150
##  7      555            600      913            854      158
##  8      557            600      709            723       53
##  9      557            600      838            846      140
## 10      558            600      753            745      138
## # ℹ 336,766 more rows
# num_range("x", 1:3)  ## x1, x2, x3にマッチする

select()で列名を変更することもできるが、明示しない変数を全て削除するのであまり便利でない。
列名の変更にはrename()を使用する。

select(flights, tail_num = tailnum)  # new_name = old_name
## # A tibble: 336,776 × 1
##    tail_num
##    <chr>   
##  1 N14228  
##  2 N24211  
##  3 N619AA  
##  4 N804JB  
##  5 N668DN  
##  6 N39463  
##  7 N516JB  
##  8 N829AS  
##  9 N593JB  
## 10 N3ALAA  
## # ℹ 336,766 more rows
rename(flights, tail_num = tailnum)
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545        -1     1004           1022
##  5  2013     1     1      554            600        -6      812            837
##  6  2013     1     1      554            558        -4      740            728
##  7  2013     1     1      555            600        -5      913            854
##  8  2013     1     1      557            600        -3      709            723
##  9  2013     1     1      557            600        -3      838            846
## 10  2013     1     1      558            600        -2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tail_num <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>

select()とヘルパー関数everything()を同時に使うことで、指定した列を先頭に移動できる。

select(flights, time_hour, air_time, everything())
## # A tibble: 336,776 × 19
##    time_hour           air_time  year month   day dep_time sched_dep_time
##    <dttm>                 <dbl> <int> <int> <int>    <int>          <int>
##  1 2013-01-01 05:00:00      227  2013     1     1      517            515
##  2 2013-01-01 05:00:00      227  2013     1     1      533            529
##  3 2013-01-01 05:00:00      160  2013     1     1      542            540
##  4 2013-01-01 05:00:00      183  2013     1     1      544            545
##  5 2013-01-01 06:00:00      116  2013     1     1      554            600
##  6 2013-01-01 05:00:00      150  2013     1     1      554            558
##  7 2013-01-01 06:00:00      158  2013     1     1      555            600
##  8 2013-01-01 06:00:00       53  2013     1     1      557            600
##  9 2013-01-01 06:00:00      140  2013     1     1      557            600
## 10 2013-01-01 06:00:00      138  2013     1     1      558            600
## # ℹ 336,766 more rows
## # ℹ 12 more variables: dep_delay <dbl>, arr_time <int>, sched_arr_time <int>,
## #   arr_delay <dbl>, carrier <chr>, flight <int>, tailnum <chr>, origin <chr>,
## #   dest <chr>, distance <dbl>, hour <dbl>, minute <dbl>

選択したい列名を文字列ベクトルで指定したい場合はone_of()を使用する。

vars <- c("year", "month", "day", "dep_delay", "arr_delay")
select(flights, one_of(vars))   ## one_of()は文字列をselectにわたす場合に必要
## # A tibble: 336,776 × 5
##     year month   day dep_delay arr_delay
##    <int> <int> <int>     <dbl>     <dbl>
##  1  2013     1     1         2        11
##  2  2013     1     1         4        20
##  3  2013     1     1         2        33
##  4  2013     1     1        -1       -18
##  5  2013     1     1        -6       -25
##  6  2013     1     1        -4        12
##  7  2013     1     1        -5        19
##  8  2013     1     1        -3       -14
##  9  2013     1     1        -3        -8
## 10  2013     1     1        -2         8
## # ℹ 336,766 more rows
# select(flights, vars)  # 直接使用すると警告が出る

3.5.4 mutate

flights_sml <- select(flights,  ## 変数が少ないデータセットを作る
  year:day, ends_with("delay"),
  distance,
  air_time)

mutate(flights_sml,
       gain = arr_delay - dep_delay,  # 新しく定義した列を最後尾に追加する
       speed = distance/air_time *60,
       hours = air_time*60,
       gain_per_hour = gain/hours) ## 作成したばかりの列を参照できる
## # A tibble: 336,776 × 11
##     year month   day dep_delay arr_delay distance air_time  gain speed hours
##    <int> <int> <int>     <dbl>     <dbl>    <dbl>    <dbl> <dbl> <dbl> <dbl>
##  1  2013     1     1         2        11     1400      227     9  370. 13620
##  2  2013     1     1         4        20     1416      227    16  374. 13620
##  3  2013     1     1         2        33     1089      160    31  408.  9600
##  4  2013     1     1        -1       -18     1576      183   -17  517. 10980
##  5  2013     1     1        -6       -25      762      116   -19  394.  6960
##  6  2013     1     1        -4        12      719      150    16  288.  9000
##  7  2013     1     1        -5        19     1065      158    24  404.  9480
##  8  2013     1     1        -3       -14      229       53   -11  259.  3180
##  9  2013     1     1        -3        -8      944      140    -5  405.  8400
## 10  2013     1     1        -2         8      733      138    10  319.  8280
## # ℹ 336,766 more rows
## # ℹ 1 more variable: gain_per_hour <dbl>

普通の書き方では新しく作る列名には変数を使用できず、変数名がそのまま新規の列名として解釈される。

flights_sml2 <- select(flights,  ## 変数が少ないデータセットを作る
  month,
  distance,
  air_time)
for (i in 1:3){
  new_col_name = paste("new_col", i, sep = "_")
  flights_sml2 <- 
    flights_sml2 %>% 
      mutate(new_col_name = n())
}
flights_sml2
## # A tibble: 336,776 × 4
##    month distance air_time new_col_name
##    <int>    <dbl>    <dbl>        <int>
##  1     1     1400      227       336776
##  2     1     1416      227       336776
##  3     1     1089      160       336776
##  4     1     1576      183       336776
##  5     1      762      116       336776
##  6     1      719      150       336776
##  7     1     1065      158       336776
##  8     1      229       53       336776
##  9     1      944      140       336776
## 10     1      733      138       336776
## # ℹ 336,766 more rows

列名に変数を使用するには先頭に!!（読み方:bang bang）を付けると変数として解釈される。またこの場合は代入に=ではなく:=を使用する必要がある。

for (i in 1:3){
  new_col_name = paste("new_col", i, sep = "_")
  flights_sml2 <- 
    flights_sml2 %>% 
      mutate(!!new_col_name := n())
}
flights_sml2
## # A tibble: 336,776 × 7
##    month distance air_time new_col_name new_col_1 new_col_2 new_col_3
##    <int>    <dbl>    <dbl>        <int>     <int>     <int>     <int>
##  1     1     1400      227       336776    336776    336776    336776
##  2     1     1416      227       336776    336776    336776    336776
##  3     1     1089      160       336776    336776    336776    336776
##  4     1     1576      183       336776    336776    336776    336776
##  5     1      762      116       336776    336776    336776    336776
##  6     1      719      150       336776    336776    336776    336776
##  7     1     1065      158       336776    336776    336776    336776
##  8     1      229       53       336776    336776    336776    336776
##  9     1      944      140       336776    336776    336776    336776
## 10     1      733      138       336776    336776    336776    336776
## # ℹ 336,766 more rows

transmute()は作成した列を抽出する

transmute(flights,
          gain = arr_delay - dep_delay,
          hours = 60*air_time,
          gain_per_hours = gain/hours) 
## # A tibble: 336,776 × 3
##     gain hours gain_per_hours
##    <dbl> <dbl>          <dbl>
##  1     9 13620       0.000661
##  2    16 13620       0.00117 
##  3    31  9600       0.00323 
##  4   -17 10980      -0.00155 
##  5   -19  6960      -0.00273 
##  6    16  9000       0.00178 
##  7    24  9480       0.00253 
##  8   -11  3180      -0.00346 
##  9    -5  8400      -0.000595
## 10    10  8280       0.00121 
## # ℹ 336,766 more rows

3.6 データの集計

3.6.1 要約

summarise()はデータフレームを要約して結果を1行で出力する。後述のgroup_by()と組み合わせて使用することが多い。
便利な要約関数

  • 中央値：median()
  • 散らばりの代表値：sd(), IOR(), mad()
  • ランクの代表値： min(), max(), quantile()
  • 位置の代表値：first(), nth(x, 2), last()
summarise(flights, mean_delay= mean(dep_delay, na.rm =T), max_delay = max(dep_delay, na.rm = T)) ## 欠損値は除く
## # A tibble: 1 × 2
##   mean_delay max_delay
##        <dbl>     <dbl>
## 1       12.6      1301

3.6.2 グループ分け

group_by()はデータを指定されたカテゴリ変数の値ごとにグループ分けし、分析単位をデータセット全体から個別グループに変更する。

flights %>% group_by(year, month) %>% 
  summarise(mean_dep_delay=mean(dep_delay, na.rm =T))  # 各年の各月ごとのdep_delayの平均値を計算
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 3
## # Groups:   year [1]
##     year month mean_dep_delay
##    <int> <int>          <dbl>
##  1  2013     1          10.0 
##  2  2013     2          10.8 
##  3  2013     3          13.2 
##  4  2013     4          13.9 
##  5  2013     5          13.0 
##  6  2013     6          20.8 
##  7  2013     7          21.7 
##  8  2013     8          12.6 
##  9  2013     9           6.72
## 10  2013    10           6.24
## 11  2013    11           5.44
## 12  2013    12          16.6
flights %>% group_by(carrier) %>% 
  filter(dep_delay >= 0) %>%  # グループごとにfilter()を適用
  summarise(mean_dep_delay=mean(dep_delay, na.rm =T))
## # A tibble: 16 × 2
##    carrier mean_dep_delay
##    <chr>            <dbl>
##  1 9E                44.9
##  2 AA                32.1
##  3 AS                27.9
##  4 B6                35.2
##  5 DL                31.5
##  6 EV                47.0
##  7 F9                40.0
##  8 FL                37.8
##  9 HA                37.3
## 10 MQ                38.6
## 11 OO                58  
## 12 UA                26.6
## 13 US                29.2
## 14 VX                29.1
## 15 WN                30.3
## 16 YV                49.2
flights %>%
  group_by(dest) %>%
  summarise(
    count = n(),
    dist = mean(distance, na.rm=T),
    delay = mean(arr_delay, na.rm=T)
    ) %>%
  filter(count >20, dest != "HML")
## # A tibble: 97 × 4
##    dest  count  dist delay
##    <chr> <int> <dbl> <dbl>
##  1 ABQ     254 1826   4.38
##  2 ACK     265  199   4.85
##  3 ALB     439  143  14.4 
##  4 ATL   17215  757. 11.3 
##  5 AUS    2439 1514.  6.02
##  6 AVL     275  584.  8.00
##  7 BDL     443  116   7.05
##  8 BGR     375  378   8.03
##  9 BHM     297  866. 16.9 
## 10 BNA    6333  758. 11.8 
## # ℹ 87 more rows

3.6.3 個数の計測

n()は各グループの行数（サイズ）を返す。非欠損値の個数を数えるにはsum(!is.na(x))を使用し、一意な値の個数を数えるにはn_distinct()を使う。

flights %>% 
  group_by(carrier, origin) %>% 
  summarise(count = n(), unique_count = n_distinct(tailnum))
## `summarise()` has grouped output by 'carrier'. You can override using the
## `.groups` argument.
## # A tibble: 35 × 4
## # Groups:   carrier [16]
##    carrier origin count unique_count
##    <chr>   <chr>  <int>        <int>
##  1 9E      EWR     1268          199
##  2 9E      JFK    14651          204
##  3 9E      LGA     2541          192
##  4 AA      EWR     3487          481
##  5 AA      JFK    13783          409
##  6 AA      LGA    15459          431
##  7 AS      EWR      714           84
##  8 B6      EWR     6557          190
##  9 B6      JFK    42076          193
## 10 B6      LGA     6002          128
## # ℹ 25 more rows

count(a, b)はdf %>% group_by(a, b) %>% summarise(count = n())に等しい。
オプションで重み変数wtを追加することもできる。

flights %>% count(year, dest)
## # A tibble: 105 × 3
##     year dest      n
##    <int> <chr> <int>
##  1  2013 ABQ     254
##  2  2013 ACK     265
##  3  2013 ALB     439
##  4  2013 ANC       8
##  5  2013 ATL   17215
##  6  2013 AUS    2439
##  7  2013 AVL     275
##  8  2013 BDL     443
##  9  2013 BGR     375
## 10  2013 BHM     297
## # ℹ 95 more rows
flights %>% count(tailnum, wt = distance)  # 個別の飛行機の総飛行距離を「カウント」できる
## # A tibble: 4,044 × 2
##    tailnum      n
##    <chr>    <dbl>
##  1 D942DN    3418
##  2 N0EGMQ  250866
##  3 N10156  115966
##  4 N102UW   25722
##  5 N103US   24619
##  6 N104UW   25157
##  7 N10575  150194
##  8 N105UW   23618
##  9 N107US   21677
## 10 N108UW   32070
## # ℹ 4,034 more rows

3.6.4 最頻値

最頻値を直接求める関数は無いので、個数や件数をカウントしてからその最大値を抽出する。
例：月ごとの目的地（dest）の最頻値を求める。
monthとdestでグループ化して個数countを求める。さらにcountの最大値max_countを求め、それとcountが一致する行を抽出。

flights %>% 
  group_by(month, dest) %>% 
  summarise(count = n()) %>% 
  mutate(max_count = max(count)) %>% 
  filter(count == max_count)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 4
## # Groups:   month [12]
##    month dest  count max_count
##    <int> <chr> <int>     <int>
##  1     1 ATL    1396      1396
##  2     2 ATL    1267      1267
##  3     3 ATL    1448      1448
##  4     4 ATL    1490      1490
##  5     5 ORD    1582      1582
##  6     6 ORD    1547      1547
##  7     7 ORD    1573      1573
##  8     8 ORD    1604      1604
##  9     9 ORD    1582      1582
## 10    10 ORD    1604      1604
## 11    11 ATL    1384      1384
## 12    12 ATL    1463      1463

countの最大値はcount %>% max()と書くこともできる。

flights %>% 
  group_by(month, dest) %>% 
  summarise(count = n()) %>% 
  filter(count == count %>% max())
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 3
## # Groups:   month [12]
##    month dest  count
##    <int> <chr> <int>
##  1     1 ATL    1396
##  2     2 ATL    1267
##  3     3 ATL    1448
##  4     4 ATL    1490
##  5     5 ORD    1582
##  6     6 ORD    1547
##  7     7 ORD    1573
##  8     8 ORD    1604
##  9     9 ORD    1582
## 10    10 ORD    1604
## 11    11 ATL    1384
## 12    12 ATL    1463

クロス集計表を使用する方法。
table()でクロス集計が可能。

table_dest <- table(flights$month, flights$dest) 
table_dest[, 1:10]
##     
##       ABQ  ACK  ALB  ANC  ATL  AUS  AVL  BDL  BGR  BHM
##   1     0    0   64    0 1396  169    2   37    0   25
##   2     0    0   58    0 1267  167    0   46    0   23
##   3     0    0   57    0 1448  273    0   62    2   26
##   4     9    0   13    0 1490  209    3   61   29   26
##   5    31   21   59    0 1499  204   35   57   29   27
##   6    30   43   34    0 1438  203   40   48   28   27
##   7    31   66   15    4 1511  213   40    2   30   29
##   8    31   67   20    4 1507  211   35   50   30   30
##   9    30   45   20    0 1364  205   30   39   56   25
##   10   31   23    1    0 1448  214   31    9   60   27
##   11   30    0   46    0 1384  186   30    3   55   22
##   12   31    0   52    0 1463  185   29   29   56   10

forループでクロス集計表のi行目の最大の要素の目的地（dest）を取り出す。which.max()は最大値のインデックスを返す。

months <- names(table_dest[,1]) # 1列目を取り出す

mode_dest <- c()
for (i in 1:length(months)){
    mode_dest[i] <- names(which.max(table_dest[i,])) # i行目のデータで最大の要素の名前を取り出す
    }
data.frame(month = months, mode_dest = mode_dest)
##    month mode_dest
## 1      1       ATL
## 2      2       ATL
## 3      3       ATL
## 4      4       ATL
## 5      5       ORD
## 6      6       ORD
## 7      7       ORD
## 8      8       ORD
## 9      9       ORD
## 10    10       ORD
## 11    11       ATL
## 12    12       ATL

3.6.5 ウインドウ関数

3.6.5.1 scoped functionとacross()

filter(), arrange(), select(), mutate(), summarise(), group_by()には複数の列を同時に処理する、scoped functionが用意されている。

  • _all: 全ての列に同じ操作を行う。
  • _if: 条件を満たす列に同じ操作を行う。
  • _at: 指定した複数の列に同じ操作を行う。
flights %>% select(ends_with("time")) %>%  # "time"で終わる列を選択
  mutate_all(abs)     # 関数abs()で全ての列の絶対値を求める
## # A tibble: 336,776 × 5
##    dep_time sched_dep_time arr_time sched_arr_time air_time
##       <int>          <int>    <int>          <int>    <dbl>
##  1      517            515      830            819      227
##  2      533            529      850            830      227
##  3      542            540      923            850      160
##  4      544            545     1004           1022      183
##  5      554            600      812            837      116
##  6      554            558      740            728      150
##  7      555            600      913            854      158
##  8      557            600      709            723       53
##  9      557            600      838            846      140
## 10      558            600      753            745      138
## # ℹ 336,766 more rows
flights %>% mutate_if(is.numeric, abs)  # 数値型のデータ列に対してのみ、abs()を適用
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545         1     1004           1022
##  5  2013     1     1      554            600         6      812            837
##  6  2013     1     1      554            558         4      740            728
##  7  2013     1     1      555            600         5      913            854
##  8  2013     1     1      557            600         3      709            723
##  9  2013     1     1      557            600         3      838            846
## 10  2013     1     1      558            600         2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% mutate_at(c("dep_delay", "arr_delay"), abs)  # 列を直接指定できる
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545         1     1004           1022
##  5  2013     1     1      554            600         6      812            837
##  6  2013     1     1      554            558         4      740            728
##  7  2013     1     1      555            600         5      913            854
##  8  2013     1     1      557            600         3      709            723
##  9  2013     1     1      557            600         3      838            846
## 10  2013     1     1      558            600         2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>

最近ではscoped functionの代わりにacross()の使用が推奨されている。

flights %>% select(ends_with("time")) %>%  # "time"で終わる列を選択
  mutate(across(everything(), abs))     # everything()で全ての列を指定
## # A tibble: 336,776 × 5
##    dep_time sched_dep_time arr_time sched_arr_time air_time
##       <int>          <int>    <int>          <int>    <dbl>
##  1      517            515      830            819      227
##  2      533            529      850            830      227
##  3      542            540      923            850      160
##  4      544            545     1004           1022      183
##  5      554            600      812            837      116
##  6      554            558      740            728      150
##  7      555            600      913            854      158
##  8      557            600      709            723       53
##  9      557            600      838            846      140
## 10      558            600      753            745      138
## # ℹ 336,766 more rows
flights %>% mutate(across(is.numeric, abs))  # 数値型のデータ列に対してのみ、abs()を適用
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(is.numeric, abs)`.
## Caused by warning:
## ! Use of bare predicate functions was deprecated in tidyselect 1.1.0.
## ℹ Please use wrap predicates in `where()` instead.
##   # Was:
##   data %>% select(is.numeric)
## 
##   # Now:
##   data %>% select(where(is.numeric))
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545         1     1004           1022
##  5  2013     1     1      554            600         6      812            837
##  6  2013     1     1      554            558         4      740            728
##  7  2013     1     1      555            600         5      913            854
##  8  2013     1     1      557            600         3      709            723
##  9  2013     1     1      557            600         3      838            846
## 10  2013     1     1      558            600         2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% mutate(across(c("dep_delay", "arr_delay"), abs))  # 列を直接指定できる
## # A tibble: 336,776 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      517            515         2      830            819
##  2  2013     1     1      533            529         4      850            830
##  3  2013     1     1      542            540         2      923            850
##  4  2013     1     1      544            545         1     1004           1022
##  5  2013     1     1      554            600         6      812            837
##  6  2013     1     1      554            558         4      740            728
##  7  2013     1     1      555            600         5      913            854
##  8  2013     1     1      557            600         3      709            723
##  9  2013     1     1      557            600         3      838            846
## 10  2013     1     1      558            600         2      753            745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% summarise(across(ends_with("time"), ~ mean(.x, na.rm = TRUE)))  # ラムダ関数も使える
## # A tibble: 1 × 5
##   dep_time sched_dep_time arr_time sched_arr_time air_time
##      <dbl>          <dbl>    <dbl>          <dbl>    <dbl>
## 1    1349.          1344.    1502.          1536.     151.

3.6.5.2 差分の計算

dplyr::lag()やdplyr::lead()により変数（ベクトル）の前後の値を取得して、差分を簡単に計算することができる。ずらすことにより欠損値になる箇所はdefaultで指定した値で埋められる。 単にlag()と書くとstats::lag()関数になるので注意。

x <- 1:10
lag1 <- dplyr::lag(x, n = 1, default = NA) # nはずらすindexの大きさ、
lead1 <- dplyr::lead(x)
data.frame(x, lag1, lead1)
##     x lag1 lead1
## 1   1   NA     2
## 2   2    1     3
## 3   3    2     4
## 4   4    3     5
## 5   5    4     6
## 6   6    5     7
## 7   7    6     8
## 8   8    7     9
## 9   9    8    10
## 10 10    9    NA

月ごとの総飛行距離の前月との差分を計算。

flights %>% 
  group_by(month) %>% 
  summarise(sum_distance = sum(distance)) %>% 
  mutate(lag_month = lag(month), lag_sum_distance = lag(sum_distance),
         diff = sum_distance - lag_sum_distance)
## # A tibble: 12 × 5
##    month sum_distance lag_month lag_sum_distance     diff
##    <int>        <dbl>     <int>            <dbl>    <dbl>
##  1     1     27188805        NA               NA       NA
##  2     2     24975509         1         27188805 -2213296
##  3     3     29179636         2         24975509  4204127
##  4     4     29427294         3         29179636   247658
##  5     5     29974128         4         29427294   546834
##  6     6     29856388         5         29974128  -117740
##  7     7     31149199         6         29856388  1292811
##  8     8     31149334         7         31149199      135
##  9     9     28711426         8         31149334 -2437908
## 10    10     30012086         9         28711426  1300660
## 11    11     28639718        10         30012086 -1372368
## 12    12     29954084        11         28639718  1314366

3.6.6 重複の削除

dplyr::distinct()によりデータフレームから重複を削除できる。引数を指定しないと重複の無いデータフレームを返す。

# originとdestの組み合わせ
flights %>% 
  select(origin, dest) %>% 
  distinct()
## # A tibble: 224 × 2
##    origin dest 
##    <chr>  <chr>
##  1 EWR    IAH  
##  2 LGA    IAH  
##  3 JFK    MIA  
##  4 JFK    BQN  
##  5 LGA    ATL  
##  6 EWR    ORD  
##  7 EWR    FLL  
##  8 LGA    IAD  
##  9 JFK    MCO  
## 10 LGA    ORD  
## # ℹ 214 more rows

引数に列名を指定することで、重複を確認する列を指定することができる。列名は複数指定することも可能。

flights %>% 
  select(origin, dest) %>% 
  distinct(origin)
## # A tibble: 3 × 1
##   origin
##   <chr> 
## 1 EWR   
## 2 LGA   
## 3 JFK
flights %>% 
  select(origin, dest) %>% 
  distinct(dest)
## # A tibble: 105 × 1
##    dest 
##    <chr>
##  1 IAH  
##  2 MIA  
##  3 BQN  
##  4 ATL  
##  5 ORD  
##  6 FLL  
##  7 IAD  
##  8 MCO  
##  9 PBI  
## 10 TPA  
## # ℹ 95 more rows
flights %>% 
  select(origin, dest) %>% 
  distinct(origin, dest) # distinct()と同じ
## # A tibble: 224 × 2
##    origin dest 
##    <chr>  <chr>
##  1 EWR    IAH  
##  2 LGA    IAH  
##  3 JFK    MIA  
##  4 JFK    BQN  
##  5 LGA    ATL  
##  6 EWR    ORD  
##  7 EWR    FLL  
##  8 LGA    IAD  
##  9 JFK    MCO  
## 10 LGA    ORD  
## # ℹ 214 more rows

引数.keep_all = TRUEとすると指定した列以外の全変数も残す。この時残るのは初出の行。

flights %>% 
  select(origin, dest) %>% 
  distinct(origin, .keep_all = TRUE)
## # A tibble: 3 × 2
##   origin dest 
##   <chr>  <chr>
## 1 EWR    IAH  
## 2 LGA    IAH  
## 3 JFK    MIA

3.7 関係データ

関係データベース（Relational database）とは表形式の複数のデータを関連付けて扱うデータベース。
表の組を結びつける変数はkeyと呼ばれる。keyは1つまたは複数で観測を一意に表し、以下の2種類がある：

  • primary key 観測を一意に識別する
  • foreign key 他の表の観測を一意に識別する
    変数がprimary keyとforeign keyの両方を兼ねる場合もある。
library(tidyverse)
library(nycflights13)
airlines
## # A tibble: 16 × 2
##    carrier name                       
##    <chr>   <chr>                      
##  1 9E      Endeavor Air Inc.          
##  2 AA      American Airlines Inc.     
##  3 AS      Alaska Airlines Inc.       
##  4 B6      JetBlue Airways            
##  5 DL      Delta Air Lines Inc.       
##  6 EV      ExpressJet Airlines Inc.   
##  7 F9      Frontier Airlines Inc.     
##  8 FL      AirTran Airways Corporation
##  9 HA      Hawaiian Airlines Inc.     
## 10 MQ      Envoy Air                  
## 11 OO      SkyWest Airlines Inc.      
## 12 UA      United Air Lines Inc.      
## 13 US      US Airways Inc.            
## 14 VX      Virgin America             
## 15 WN      Southwest Airlines Co.     
## 16 YV      Mesa Airlines Inc.
airports
## # A tibble: 1,458 × 8
##    faa   name                             lat    lon   alt    tz dst   tzone    
##    <chr> <chr>                          <dbl>  <dbl> <dbl> <dbl> <chr> <chr>    
##  1 04G   Lansdowne Airport               41.1  -80.6  1044    -5 A     America/…
##  2 06A   Moton Field Municipal Airport   32.5  -85.7   264    -6 A     America/…
##  3 06C   Schaumburg Regional             42.0  -88.1   801    -6 A     America/…
##  4 06N   Randall Airport                 41.4  -74.4   523    -5 A     America/…
##  5 09J   Jekyll Island Airport           31.1  -81.4    11    -5 A     America/…
##  6 0A9   Elizabethton Municipal Airport  36.4  -82.2  1593    -5 A     America/…
##  7 0G6   Williams County Airport         41.5  -84.5   730    -5 A     America/…
##  8 0G7   Finger Lakes Regional Airport   42.9  -76.8   492    -5 A     America/…
##  9 0P2   Shoestring Aviation Airfield    39.8  -76.6  1000    -5 U     America/…
## 10 0S9   Jefferson County Intl           48.1 -123.    108    -8 A     America/…
## # ℹ 1,448 more rows
planes
## # A tibble: 3,322 × 9
##    tailnum  year type              manufacturer model engines seats speed engine
##    <chr>   <int> <chr>             <chr>        <chr>   <int> <int> <int> <chr> 
##  1 N10156   2004 Fixed wing multi… EMBRAER      EMB-…       2    55    NA Turbo…
##  2 N102UW   1998 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  3 N103US   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  4 N104UW   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  5 N10575   2002 Fixed wing multi… EMBRAER      EMB-…       2    55    NA Turbo…
##  6 N105UW   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  7 N107US   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  8 N108UW   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
##  9 N109UW   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
## 10 N110UW   1999 Fixed wing multi… AIRBUS INDU… A320…       2   182    NA Turbo…
## # ℹ 3,312 more rows
weather
## # A tibble: 26,115 × 15
##    origin  year month   day  hour  temp  dewp humid wind_dir wind_speed
##    <chr>  <int> <int> <int> <int> <dbl> <dbl> <dbl>    <dbl>      <dbl>
##  1 EWR     2013     1     1     1  39.0  26.1  59.4      270      10.4 
##  2 EWR     2013     1     1     2  39.0  27.0  61.6      250       8.06
##  3 EWR     2013     1     1     3  39.0  28.0  64.4      240      11.5 
##  4 EWR     2013     1     1     4  39.9  28.0  62.2      250      12.7 
##  5 EWR     2013     1     1     5  39.0  28.0  64.4      260      12.7 
##  6 EWR     2013     1     1     6  37.9  28.0  67.2      240      11.5 
##  7 EWR     2013     1     1     7  39.0  28.0  64.4      240      15.0 
##  8 EWR     2013     1     1     8  39.9  28.0  62.2      250      10.4 
##  9 EWR     2013     1     1     9  39.9  28.0  62.2      260      15.0 
## 10 EWR     2013     1     1    10  41    28.0  59.6      260      13.8 
## # ℹ 26,105 more rows
## # ℹ 5 more variables: wind_gust <dbl>, precip <dbl>, pressure <dbl>,
## #   visib <dbl>, time_hour <dttm>

primary keyが本当に観測を一意に識別しているか検証する。primary keyが無い場合はmutate()とrow_number()等で代替keyを作る。

planes %>% count(tailnum) %>% filter(n > 1)  # primary keyのcount()で、各エントリについてnが1より大きくないか調べる
## # A tibble: 0 × 2
## # ℹ 2 variables: tailnum <chr>, n <int>
# weatherの観測を一意に識別するには、5つの変数が必要, countは自動的にgroup_by()することに注意
weather %>% count(year, month, day, hour, origin) %>% filter(n > 1)
## # A tibble: 3 × 6
##    year month   day  hour origin     n
##   <int> <int> <int> <int> <chr>  <int>
## 1  2013    11     3     1 EWR        2
## 2  2013    11     3     1 JFK        2
## 3  2013    11     3     1 LGA        2

3.7.1 更新ジョイン

更新ジョイン mutating joinは2つの表（X, Y）の変数を組み合わせる。 keyとマッチする観測を探し出し、変数を１つの表からもう１つの表にコピーする。mutate()と同様、ジョイン関数は変数を右側に追加する。

flights2 <- flights %>% select(year:day, hour, origin, dest, tailnum, carrier)
flights2
## # A tibble: 336,776 × 8
##     year month   day  hour origin dest  tailnum carrier
##    <int> <int> <int> <dbl> <chr>  <chr> <chr>   <chr>  
##  1  2013     1     1     5 EWR    IAH   N14228  UA     
##  2  2013     1     1     5 LGA    IAH   N24211  UA     
##  3  2013     1     1     5 JFK    MIA   N619AA  AA     
##  4  2013     1     1     5 JFK    BQN   N804JB  B6     
##  5  2013     1     1     6 LGA    ATL   N668DN  DL     
##  6  2013     1     1     5 EWR    ORD   N39463  UA     
##  7  2013     1     1     6 EWR    FLL   N516JB  B6     
##  8  2013     1     1     6 LGA    IAD   N829AS  EV     
##  9  2013     1     1     6 JFK    MCO   N593JB  B6     
## 10  2013     1     1     6 LGA    ORD   N3ALAA  AA     
## # ℹ 336,766 more rows
# flightsに航空会社の正式名を追加するため、airlinesとflightsをleft_joinで組み合わせる
flights2 %>% select(-origin, -dest) %>% left_join(airlines, by = "carrier")
## # A tibble: 336,776 × 7
##     year month   day  hour tailnum carrier name                    
##    <int> <int> <int> <dbl> <chr>   <chr>   <chr>                   
##  1  2013     1     1     5 N14228  UA      United Air Lines Inc.   
##  2  2013     1     1     5 N24211  UA      United Air Lines Inc.   
##  3  2013     1     1     5 N619AA  AA      American Airlines Inc.  
##  4  2013     1     1     5 N804JB  B6      JetBlue Airways         
##  5  2013     1     1     6 N668DN  DL      Delta Air Lines Inc.    
##  6  2013     1     1     5 N39463  UA      United Air Lines Inc.   
##  7  2013     1     1     6 N516JB  B6      JetBlue Airways         
##  8  2013     1     1     6 N829AS  EV      ExpressJet Airlines Inc.
##  9  2013     1     1     6 N593JB  B6      JetBlue Airways         
## 10  2013     1     1     6 N3ALAA  AA      American Airlines Inc.  
## # ℹ 336,766 more rows

3.7.1.1 内部ジョイン

内部ジョイン inner join, keyが共通した行が抽出される、 X AND Y。一般に観測が失われるため分析には不向き。

x <- tribble(
  ~key, ~val_x,
  1, "x1",
  2, "x2",
  3, "x3"
)
y <- tribble(
  ~key, ~val_y,
  1, "y1",
  2, "y2",
  4, "y3"
)

x %>% 
  inner_join(y, by = "key")
## # A tibble: 2 × 3
##     key val_x val_y
##   <dbl> <chr> <chr>
## 1     1 x1    y1   
## 2     2 x2    y2

3.7.1.2 外部ジョイン

内部ジョインが両方の表にある観測を保持する一方、外部ジョイン outer joinは、少なくとも1つの表にある観測を保持する。以下の3種類がある：

  1. left join：Xの全観測を保持する、マッチが無くても元の観測が保持されるので一番良く使われる
  2. right join：Yの全観測を保持する
  3. full join：XとY全ての観測を保持する、X OR Y

重複キー

# 1つの表に重複キーがある場合
x <- tribble(
  ~key, ~val_x,
  1, "x1",
  2, "x2",
  2, "x3",
  1, "x4"
)
y <- tribble(
  ~key, ~val_y,
  1, "y1",
  2, "y2"
)
left_join(x, y, by = "key")
## # A tibble: 4 × 3
##     key val_x val_y
##   <dbl> <chr> <chr>
## 1     1 x1    y1   
## 2     2 x2    y2   
## 3     2 x3    y2   
## 4     1 x4    y1
# 両方の表に重複キーが含まれる場合、どちらの表でも観測を一意に識別できないので、普通はエラーになる。joinすると全ての組み合わせが得られる
x <- tribble(
  ~key, ~val_x,
  1, "x1",
  2, "x2",
  2, "x3",
  3, "x4"
)
y <- tribble(
  ~key, ~val_y,
  1, "y1",
  2, "y2",
  2, "y3",
  3, "y4"
)
left_join(x, y, by = "key")
## Warning in left_join(x, y, by = "key"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 2 of `x` matches multiple rows in `y`.
## ℹ Row 2 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## # A tibble: 6 × 3
##     key val_x val_y
##   <dbl> <chr> <chr>
## 1     1 x1    y1   
## 2     2 x2    y2   
## 3     2 x2    y3   
## 4     2 x3    y2   
## 5     2 x3    y3   
## 6     3 x4    y4

Key列の指定
defaultではby =NULLで、両方の表に現れる全ての変数を使い、natural joinする

flights2 %>% 
  left_join(weather)
## Joining with `by = join_by(year, month, day, hour, origin)`
## # A tibble: 336,776 × 18
##     year month   day  hour origin dest  tailnum carrier  temp  dewp humid
##    <int> <int> <int> <dbl> <chr>  <chr> <chr>   <chr>   <dbl> <dbl> <dbl>
##  1  2013     1     1     5 EWR    IAH   N14228  UA       39.0  28.0  64.4
##  2  2013     1     1     5 LGA    IAH   N24211  UA       39.9  25.0  54.8
##  3  2013     1     1     5 JFK    MIA   N619AA  AA       39.0  27.0  61.6
##  4  2013     1     1     5 JFK    BQN   N804JB  B6       39.0  27.0  61.6
##  5  2013     1     1     6 LGA    ATL   N668DN  DL       39.9  25.0  54.8
##  6  2013     1     1     5 EWR    ORD   N39463  UA       39.0  28.0  64.4
##  7  2013     1     1     6 EWR    FLL   N516JB  B6       37.9  28.0  67.2
##  8  2013     1     1     6 LGA    IAD   N829AS  EV       39.9  25.0  54.8
##  9  2013     1     1     6 JFK    MCO   N593JB  B6       37.9  27.0  64.3
## 10  2013     1     1     6 LGA    ORD   N3ALAA  AA       39.9  25.0  54.8
## # ℹ 336,766 more rows
## # ℹ 7 more variables: wind_dir <dbl>, wind_speed <dbl>, wind_gust <dbl>,
## #   precip <dbl>, pressure <dbl>, visib <dbl>, time_hour <dttm>
# flight2とplanesにはどちらもyearがあるが、それぞれ異なる意味なのでtailnumだけでjoinしたい

flights2 %>% 
  left_join(planes, by = "tailnum")
## # A tibble: 336,776 × 16
##    year.x month   day  hour origin dest  tailnum carrier year.y type            
##     <int> <int> <int> <dbl> <chr>  <chr> <chr>   <chr>    <int> <chr>           
##  1   2013     1     1     5 EWR    IAH   N14228  UA        1999 Fixed wing mult…
##  2   2013     1     1     5 LGA    IAH   N24211  UA        1998 Fixed wing mult…
##  3   2013     1     1     5 JFK    MIA   N619AA  AA        1990 Fixed wing mult…
##  4   2013     1     1     5 JFK    BQN   N804JB  B6        2012 Fixed wing mult…
##  5   2013     1     1     6 LGA    ATL   N668DN  DL        1991 Fixed wing mult…
##  6   2013     1     1     5 EWR    ORD   N39463  UA        2012 Fixed wing mult…
##  7   2013     1     1     6 EWR    FLL   N516JB  B6        2000 Fixed wing mult…
##  8   2013     1     1     6 LGA    IAD   N829AS  EV        1998 Fixed wing mult…
##  9   2013     1     1     6 JFK    MCO   N593JB  B6        2004 Fixed wing mult…
## 10   2013     1     1     6 LGA    ORD   N3ALAA  AA          NA <NA>            
## # ℹ 336,766 more rows
## # ℹ 6 more variables: manufacturer <chr>, model <chr>, engines <int>,
## #   seats <int>, speed <int>, engine <chr>
# この場合、出力ではyear.x、year.yと区別される

# c("a" = "b")とすると、表Xのaと表Yのbをマッチさせる
flights2 %>% 
  left_join(airports, c("dest" = "faa"))
## # A tibble: 336,776 × 15
##     year month   day  hour origin dest  tailnum carrier name     lat   lon   alt
##    <int> <int> <int> <dbl> <chr>  <chr> <chr>   <chr>   <chr>  <dbl> <dbl> <dbl>
##  1  2013     1     1     5 EWR    IAH   N14228  UA      Georg…  30.0 -95.3    97
##  2  2013     1     1     5 LGA    IAH   N24211  UA      Georg…  30.0 -95.3    97
##  3  2013     1     1     5 JFK    MIA   N619AA  AA      Miami…  25.8 -80.3     8
##  4  2013     1     1     5 JFK    BQN   N804JB  B6      <NA>    NA    NA      NA
##  5  2013     1     1     6 LGA    ATL   N668DN  DL      Harts…  33.6 -84.4  1026
##  6  2013     1     1     5 EWR    ORD   N39463  UA      Chica…  42.0 -87.9   668
##  7  2013     1     1     6 EWR    FLL   N516JB  B6      Fort …  26.1 -80.2     9
##  8  2013     1     1     6 LGA    IAD   N829AS  EV      Washi…  38.9 -77.5   313
##  9  2013     1     1     6 JFK    MCO   N593JB  B6      Orlan…  28.4 -81.3    96
## 10  2013     1     1     6 LGA    ORD   N3ALAA  AA      Chica…  42.0 -87.9   668
## # ℹ 336,766 more rows
## # ℹ 3 more variables: tz <dbl>, dst <chr>, tzone <chr>
flights2 %>% 
  left_join(airports, c("origin" = "faa"))
## # A tibble: 336,776 × 15
##     year month   day  hour origin dest  tailnum carrier name     lat   lon   alt
##    <int> <int> <int> <dbl> <chr>  <chr> <chr>   <chr>   <chr>  <dbl> <dbl> <dbl>
##  1  2013     1     1     5 EWR    IAH   N14228  UA      Newar…  40.7 -74.2    18
##  2  2013     1     1     5 LGA    IAH   N24211  UA      La Gu…  40.8 -73.9    22
##  3  2013     1     1     5 JFK    MIA   N619AA  AA      John …  40.6 -73.8    13
##  4  2013     1     1     5 JFK    BQN   N804JB  B6      John …  40.6 -73.8    13
##  5  2013     1     1     6 LGA    ATL   N668DN  DL      La Gu…  40.8 -73.9    22
##  6  2013     1     1     5 EWR    ORD   N39463  UA      Newar…  40.7 -74.2    18
##  7  2013     1     1     6 EWR    FLL   N516JB  B6      Newar…  40.7 -74.2    18
##  8  2013     1     1     6 LGA    IAD   N829AS  EV      La Gu…  40.8 -73.9    22
##  9  2013     1     1     6 JFK    MCO   N593JB  B6      John …  40.6 -73.8    13
## 10  2013     1     1     6 LGA    ORD   N3ALAA  AA      La Gu…  40.8 -73.9    22
## # ℹ 336,766 more rows
## # ℹ 3 more variables: tz <dbl>, dst <chr>, tzone <chr>

3.7.2 フィルタジョイン

Filtering join フィルタジョイン、変数ではなく観測に影響を及ぼす。

  • semi_join(x, y):セミジョインはyのとマッチするxの全ての観測を保持する。
  • anti_join(x, y):アンチジョインはyのとマッチするxの全ての観測を取り除く。
# 例） 人々が最もよく行く旅行先のトップ10を探す。
top_dest <- flights %>% count(dest, sort = T) %>% head(10)
top_dest
## # A tibble: 10 × 2
##    dest      n
##    <chr> <int>
##  1 ORD   17283
##  2 ATL   17215
##  3 LAX   16174
##  4 BOS   15508
##  5 MCO   14082
##  6 CLT   14064
##  7 SFO   13331
##  8 FLL   12055
##  9 MIA   11728
## 10 DCA    9705
# これらへのフライトを探す
flights %>% filter(dest %in% top_dest$dest)
## # A tibble: 141,145 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      542            540         2      923            850
##  2  2013     1     1      554            600        -6      812            837
##  3  2013     1     1      554            558        -4      740            728
##  4  2013     1     1      555            600        -5      913            854
##  5  2013     1     1      557            600        -3      838            846
##  6  2013     1     1      558            600        -2      753            745
##  7  2013     1     1      558            600        -2      924            917
##  8  2013     1     1      558            600        -2      923            937
##  9  2013     1     1      559            559         0      702            706
## 10  2013     1     1      600            600         0      851            858
## # ℹ 141,135 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>
# しかしこの方式を複数の変数について拡張するのは容易ではないので、代わりにsemi_joinを使用する
flights %>% semi_join(top_dest)
## Joining with `by = join_by(dest)`
## # A tibble: 141,145 × 19
##     year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
##    <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
##  1  2013     1     1      542            540         2      923            850
##  2  2013     1     1      554            600        -6      812            837
##  3  2013     1     1      554            558        -4      740            728
##  4  2013     1     1      555            600        -5      913            854
##  5  2013     1     1      557            600        -3      838            846
##  6  2013     1     1      558            600        -2      753            745
##  7  2013     1     1      558            600        -2      924            917
##  8  2013     1     1      558            600        -2      923            937
##  9  2013     1     1      559            559         0      702            706
## 10  2013     1     1      600            600         0      851            858
## # ℹ 141,135 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## #   tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## #   hour <dbl>, minute <dbl>, time_hour <dttm>

アンチジョインはジョインのミスマッチの診断に使うことができる。

# flightsとplanesを連結する時に、planesとマッチがとれないflightsが多いかどうか調べる
flights %>% anti_join(planes, by = "tailnum") %>% count(tailnum, sort = T)
## # A tibble: 722 × 2
##    tailnum     n
##    <chr>   <int>
##  1 <NA>     2512
##  2 N725MQ    575
##  3 N722MQ    513
##  4 N723MQ    507
##  5 N713MQ    483
##  6 N735MQ    396
##  7 N0EGMQ    371
##  8 N534MQ    364
##  9 N542MQ    363
## 10 N531MQ    349
## # ℹ 712 more rows

3.7.3 ジョインの問題

  1. 各表の主キーとなる変数を確認する
  2. 主キーの変数に欠損値が無いことを確認する
  3. 外部キーが他の表の主キーとマッチすることを確かめる

3.7.4 集合演算

intersect(x,y)：xとyに共通な観測だけを返す。
union(x,y)：xとyにある観測を一意にして返す。
setdiff(x,y)：xにはあるがyにはない観測を返す。

df1 <- tribble(
  ~x, ~y,
  1,  1,
  2,  1
)
df2 <- tribble(
  ~x, ~y,
  1,  1,
  1,  2
)

intersect(df1, df2)
## # A tibble: 1 × 2
##       x     y
##   <dbl> <dbl>
## 1     1     1
union(df1, df2)
## # A tibble: 3 × 2
##       x     y
##   <dbl> <dbl>
## 1     1     1
## 2     2     1
## 3     1     2
setdiff(df1, df2)
## # A tibble: 1 × 2
##       x     y
##   <dbl> <dbl>
## 1     2     1
setdiff(df2, df1)
## # A tibble: 1 × 2
##       x     y
##   <dbl> <dbl>
## 1     1     2

3.8 データの畳み込み

tibbleにはデータフレームやtibbleも格納できる。tidyr::nest()を使うと、データをtibble内に畳み込むことができる。畳み込まれたサブデータはデフォルトではリスト列”data”に格納される。名前を変更する場合は引数.keyで指定する。畳み込まれたデータ構造のことをネストデータあるいは入れ子データと呼ぶ。

# ダミーデータの生成, LETTERSとlettersは共にbuilt-in constantsでそれぞれ大文字と小文字のアルファベットが格納されている
set.seed(71)
N <- 15
dat <- tibble(tag1 = sample(LETTERS[1:3], N, replace = TRUE),
              tag2 = sample(letters[1:5], N, replace = TRUE),
              y = rnorm(N),
              x = runif(N))

tibbleを要素に持つ、リスト形式で格納されている。

dat %>% 
  nest(.key = "data")
## # A tibble: 1 × 1
##   data             
##   <list>           
## 1 <tibble [15 × 4]>

事前に変数tag1でグループ化しておくと、tag1の水準ごとにデータを畳み込む。

dat %>% 
  group_by(tag1) %>% 
  nest()
## # A tibble: 3 × 2
## # Groups:   tag1 [3]
##   tag1  data            
##   <chr> <list>          
## 1 C     <tibble [6 × 3]>
## 2 B     <tibble [5 × 3]>
## 3 A     <tibble [4 × 3]>

group_nest()を使うとグループ化と畳み込みを同時に行う。

dat %>% 
  group_nest(tag1, tag2)
## # A tibble: 9 × 3
##   tag1  tag2                data
##   <chr> <chr> <list<tibble[,2]>>
## 1 A     a                [1 × 2]
## 2 A     c                [1 × 2]
## 3 A     d                [2 × 2]
## 4 B     b                [1 × 2]
## 5 B     c                [2 × 2]
## 6 B     e                [2 × 2]
## 7 C     b                [2 × 2]
## 8 C     d                [1 × 2]
## 9 C     e                [3 × 2]

畳み込まれたデータを広げるにはunnest()を使う。リスト列が複数ある場合、広げるリスト列を()内で指定できる。

dat %>% 
  group_nest(tag1, tag2) %>% 
  unnest()
## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data)`.
## # A tibble: 15 × 4
##    tag1  tag2        y      x
##    <chr> <chr>   <dbl>  <dbl>
##  1 A     a      0.324  0.0138
##  2 A     c     -0.904  0.414 
##  3 A     d      0.900  0.427 
##  4 A     d      1.29   0.555 
##  5 B     b     -1.62   0.473 
##  6 B     c     -0.730  0.161 
##  7 B     c      1.70   0.889 
##  8 B     e      2.37   0.239 
##  9 B     e     -0.673  0.648 
## 10 C     b      0.252  0.550 
## 11 C     b      0.324  0.397 
## 12 C     d     -1.15   0.449 
## 13 C     e      0.0100 0.808 
## 14 C     e     -0.887  0.0279
## 15 C     e     -1.79   0.110

このような入れ子データ構造は、グループ分けされたサブデータに対してmap関数で統計処理を適用するのに有用。
data列の各データのxの平均値を計算して新しい列に格納してみる。map()だとリストが返り値になるが、map_dbl()だと実数が返り値。

dat_nest <- dat %>% 
  group_nest(tag1)
dat_nest %>% 
  mutate(x_mean = map_dbl(data, ~ mean(.$x)))  # .でリスト列中の各リストを表すためには~による記法が必要
## # A tibble: 3 × 3
##   tag1                data x_mean
##   <chr> <list<tibble[,3]>>  <dbl>
## 1 A                [4 × 3]  0.352
## 2 B                [5 × 3]  0.482
## 3 C                [6 × 3]  0.390

group_by()とsummarise()を使用した場合と比較して、統計要約量を別のデータフレームに格納する必要が無いのが利点。同様に統計検定も可能。

dat_nest %>%
  mutate(ttest = map(data, ~ t.test(.$x, .$y)),
         pval = map_dbl(ttest, ~.$p.value))
## # A tibble: 3 × 4
##   tag1                data ttest     pval
##   <chr> <list<tibble[,3]>> <list>   <dbl>
## 1 A                [4 × 3] <htest> 0.926 
## 2 B                [5 × 3] <htest> 0.745 
## 3 C                [6 × 3] <htest> 0.0456

各サブデータのさらに一部を抽出して検定することも可能。以下ではtag2 == "b"でないデータについて検定している。

dat_nest %>%
  mutate(x = map(data, ~ filter(., tag2 != "b") %>% .$x),
         y = map(data, ~ filter(., tag2 != "b") %>% .$y),
         ttest = map2(x, y, ~ t.test(.x, .y)),
         pval = map_dbl(ttest, ~ .$p.value))
## # A tibble: 3 × 6
##   tag1                data x         y         ttest     pval
##   <chr> <list<tibble[,3]>> <list>    <list>    <list>   <dbl>
## 1 A                [4 × 3] <dbl [4]> <dbl [4]> <htest> 0.926 
## 2 B                [5 × 3] <dbl [4]> <dbl [4]> <htest> 0.837 
## 3 C                [6 × 3] <dbl [4]> <dbl [4]> <htest> 0.0311

4 データの可視化

4.1 ggplot2によるグラフ作成

4.1.1 基本的な考え方

ggplot2はHadely Wickhamによって開発されたデータ可視化のためのパッケージ。「Grammer of Graphics」の思想に基づいている。グラフの構成要素のそれぞれを1つの層（レイヤー）と捉え、これらを重ねることでグラフを作成する。

4.1.2 用語の定義

  • データ : 視覚化される情報。データは変数とその値で構成され、データフレームにロング形式で格納される。カテゴリ変数の値は離散的であり、数値変数の値は連続的。
  • 幾何オブジェクト（geometric object） : データを表現するための散布図、ヒストグラム等のグラフ。折れ線グラフならgeom_line()といったように、幾何オブジェクトはgeom_で始まる関数で表される。
  • エステティック属性（aesthetic） : 幾何オブジェクトの視覚的性質（x, yの位置、線の色、点の形など）。エステティック属性は幾何オブジェクトごとに異なる。
  • マッピング : データの変数とエステティック属性を対応させること。mapping = aes(マッピング)で指定する。mapping=は省略可能。
  • 設定 : エステティック属性を、データ変数の値とは無関係な、ある一定の値に設定すること。
  • スケール : データ空間の値からエステティック空間へ、値をマッピングするやり方をコントロールする。
  • ガイド : 視覚属性がデータ空間にどのようにマッピングし直されるかを示す。例: 軸上の目盛りとラベル、凡例など。

4.1.3 グラフ作成の文法

まずデータの入ったデータフレームを呼び出し、ggplot()関数でデフォルトのマッピングを指定する。下の例ではxに変数coutryを、yに変数casesを指定。次にどのような幾何オブジェクトでプロットするかを指定する。下の例ではgeom_point()で散布図にする。またgeom_point()内にaes()を置き、変数yearを点の色であるcolorにマッピングする。このcolor = yearというマッピングはgeom_point()に限定されており、他の幾何オブジェクトには作用しない。もしggplot()内のaes()でcolor = yearとマッピングすれば、以下の全ての幾何オブジェクトに作用する。

library(tidyverse)
tb1<- table1 # tidyverseに収録されているデータ
plot1 <- tb1 %>% 
  ggplot(mapping = aes(x = country, y = cases)) + 
  geom_point(aes(color = year))
plot1

マッピングと設定の違い。colorの値をaes()の外でredに指定する。データの値とは無関係に全ての点が赤色になる。

plot2 <- tb1 %>% 
  ggplot(mapping = aes(x = country, y = cases)) + 
  geom_point(color = "red")
plot2

デフォルトでは変数yearは数値型であり連続的変数としてcolorにマッピングされた。factor型に変換してから離散的変数としてマッピングすることもできる。

plot3 <- tb1 %>% 
  mutate(year = as.factor(year)) %>% 
  ggplot(mapping = aes(x = country, y = cases)) + 
  geom_point(aes(color = year))
plot3

スケール変更の例。scale_color_manual()で色のスケール（色の対応のさせ方）を変更。またscale_y_continuous()でy軸の範囲を変更する。

plot3r <- tb1 %>% 
  mutate(year = as.factor(year)) %>% 
  ggplot(mapping = aes(x = country, y = cases)) + 
  geom_point(aes(color = year)) +
  scale_colour_manual(values = c("orange", "forestgreen"))+
  scale_y_continuous(limits = c(0, 250000))
  
plot3r

4.2 Rにおける色の使い方

ggplot2において色はエステティック属性の1つであり、x座標やサイズ（size）と同じように扱われる。
幾何オブジェクトの色を設定するには、geom_**()関数のaes()の 外でcolorまたはfillの値を設定する。
変数を幾何オブジェクトの色にマッピングするには、aes()の中でcolorまたはfillの値に変数をマッピングする。

library(RColorBrewer)
display.brewer.all(type = "div")

display.brewer.all(type = "seq")

display.brewer.all(type = "qual")

display.brewer.all(colorblindFriendly = TRUE)

brewer.pal(4, "Set2")
## [1] "#66C2A5" "#FC8D62" "#8DA0CB" "#E78AC3"
dput(brewer.pal(4, "Set2")) # dput()はオブジェクトをテキスト形式でコンソールに出力する、ベクトルを生成するのに便利
## c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3")
display.brewer.pal(4, "Set2")

brewer.pal.info
##          maxcolors category colorblind
## BrBG            11      div       TRUE
## PiYG            11      div       TRUE
## PRGn            11      div       TRUE
## PuOr            11      div       TRUE
## RdBu            11      div       TRUE
## RdGy            11      div      FALSE
## RdYlBu          11      div       TRUE
## RdYlGn          11      div      FALSE
## Spectral        11      div      FALSE
## Accent           8     qual      FALSE
## Dark2            8     qual       TRUE
## Paired          12     qual       TRUE
## Pastel1          9     qual      FALSE
## Pastel2          8     qual      FALSE
## Set1             9     qual      FALSE
## Set2             8     qual       TRUE
## Set3            12     qual      FALSE
## Blues            9      seq       TRUE
## BuGn             9      seq       TRUE
## BuPu             9      seq       TRUE
## GnBu             9      seq       TRUE
## Greens           9      seq       TRUE
## Greys            9      seq       TRUE
## Oranges          9      seq       TRUE
## OrRd             9      seq       TRUE
## PuBu             9      seq       TRUE
## PuBuGn           9      seq       TRUE
## PuRd             9      seq       TRUE
## Purples          9      seq       TRUE
## RdPu             9      seq       TRUE
## Reds             9      seq       TRUE
## YlGn             9      seq       TRUE
## YlGnBu           9      seq       TRUE
## YlOrBr           9      seq       TRUE
## YlOrRd           9      seq       TRUE

ggplot2の標準のカラーパレット

dum_dat <- data.frame(key=c("A", "B", "C", "D", "E", "F", "G", "H"),
                             value = c(1, 1, 1, 1, 1, 1, 1, 1), stringsAsFactors = TRUE)
std_col_pal <- dum_dat %>% 
  ggplot(aes(x = key, y = value, fill=key))+
  geom_col()+
  scale_fill_manual(values = c("#F8766D", "#CD9600", "#7CAE00", "#00BE67", "#00BFC4", "#00A9FF", "#C77CFF", "#FF61CC"))+
  xlab("")+
  ylab("")
std_col_pal

無難な3色

acceptable_3col_pal <- dum_dat %>%
  filter(key %in% c("A", "B", "C")) %>% 
  ggplot(aes(x = key, y = value, fill=key))+
  geom_col()+
  scale_fill_manual(values = c('#bdbdbd','#3182bd',"#e6550d"))+
  xlab("")+
  ylab("")
acceptable_3col_pal

無難な3色（やや明るめ）

acceptable_3col_pal <- dum_dat %>%
  filter(key %in% c("A", "B", "C")) %>% 
  ggplot(aes(x = key, y = value, fill=key))+
  geom_col()+
  # scale_fill_manual(values = c("#c1afa8", "#7ec8ef","#F8766D"))+
    # scale_fill_manual(values = c("#a9a9ac","#5da6f4", "#ff8078"))+
     # scale_fill_manual(values = c("#c1afa8","#5da6f4", "#ff8078"))+
   scale_fill_manual(values = c("#bdbdbd","#5da6f4", "#F8766D"))+
  xlab("")+
  ylab("")
acceptable_3col_pal

無難な5色

acceptable_5col_pal <- dum_dat %>%
  filter(key %in% c("A", "B", "C", "D", "E")) %>% 
  ggplot(aes(x = key, y = value, fill=key))+
  geom_col()+
  scale_fill_manual(values = c('#bdbdbd',"#F8766D", "#00BFC4", "#619CFF", "#F564E3"))+
  xlab("")+
  ylab("")
acceptable_5col_pal

5 モデル

5.1 モデル化の基本

モデルの目的は、データセットの簡単な低次元の要約。仮説生成と仮説確認には別のデータセットを使用する。
1. モデルのfamilyを定義する
2. 適合モデル（a fitted model）をデータに最も近いファミリーから探索し生成する、適合モデルは何らかの基準で「最良」のモデル

パッケージmodelrは基本Rのモデル化機能をパイプで自然に使えるようwrapしたもの。

library(tidyverse)
library(modelr) 
options(na.action = na.warn)

sim1はmodelrに格納されている模擬データセットの1つ、線形パターンy = a_0 + a_1 * xを示す。

ggplot(sim1, aes(x,y)) +
  geom_point()

sim1
## # A tibble: 30 × 2
##        x     y
##    <int> <dbl>
##  1     1  4.20
##  2     1  7.51
##  3     1  2.13
##  4     2  8.99
##  5     2 10.2 
##  6     2 11.3 
##  7     3  7.36
##  8     3 10.5 
##  9     3 10.5 
## 10     4 12.4 
## # ℹ 20 more rows

ランダムに線を重ねる。

models <- tibble(
  # -20〜40の間で250個の実数をランダムに生成
  a1 = runif(250, -20, 40),
  a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) + 
  # 250組の(a1, a2)により250本の直線が引かれる
  geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +  
  geom_point()

距離計算
model1はモデルのファミリーを表す、モデルの引数（xの値）とデータを入力にとり、モデルの予測値を出力。

model1 <- function(a, data){
  a[1] + data$x * a[2]  ##data$xはデータフレームdataの列xを表す
}
model1(c(7, 1.5), sim1)
##  [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5
## [16] 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0 22.0 22.0

モデルによるy値（予測）と実際のデータのy値（応答）の差 = 距離を計算する。複数の個別の距離を全体的な距離にまとめる必要がある、これには「偏差の二乗平均平方根 root-mean squared deviation」（RMSD）を使う。実際と予測の差を計算して、二乗し、平均を取り、その平方根を求める。

measure_distance <- function(mod, data){
  diff <- data$y - model1(mod, data)
  sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
## [1] 2.665212

modelsの中の250組全ての(a1, a2)に対してRMSDを計算。

sim1_dist <- function(a1, a2){
  measure_distance(c(a1, a2), sim1)
}
# map2は引数を2つ取ることができる、計算したRMSDを列distに格納
models <- models %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist)) 
models
## # A tibble: 250 × 3
##        a1     a2  dist
##     <dbl>  <dbl> <dbl>
##  1 -10.6  -1.67  36.9 
##  2 -15.4  -0.839 36.5 
##  3   9.19 -2.13  21.8 
##  4  -3.37  4.42   8.94
##  5 -13.6  -1.86  40.9 
##  6 -14.6   2.27  17.8 
##  7  38.9  -4.28  18.3 
##  8  30.1   1.15  21.1 
##  9  33.5   0.688 22.2 
## 10   1.85 -2.52  30.6 
## # ℹ 240 more rows

最良の10モデルをプロット、distが小さいほど良いモデル。

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), ## color = -distで最小のdistが一番明るい色になる
    data = filter(models, rank(dist) <= 10)  ## distの小さいもの10位を選ぶ
  )

モデルを観測と考えて、散布図を描き, -distで色付け、最良の10モデルを赤丸で囲む。

ggplot(models, aes(a1, a2)) +
  geom_point(
    data = filter(models, rank(dist) <= 10),
    size =4, color = "red"
  ) +
  geom_point(aes(color = -dist))

格子探索、より系統的に等間隔に格子点を生成する。どこに最良のモデルがあるか見当がついたので、大雑把に格子のパラメータを選ぶ。expand.grid関数はn個のベクトルに含まれる要素の全ての組み合わせをすばやく書き出す。

grid <- expand.grid(
  a1 = seq(-5, 20, length = 25),
  a2 = seq(1, 3, length = 25)
) %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))

grid %>% 
  ggplot(aes(a1, a2)) +
  geom_point(
    data = filter(grid, rank(dist) <= 10),
    size = 4, color = "red"
  ) +
  geom_point(aes(color = -dist))

最良の10モデルを元のデータに重ねる。

ggplot(sim1, aes(x, y)) +
  geom_point(size = 2, color = "grey30") +
  geom_abline(
    aes(intercept = a1, slope = a2, color = -dist),
    data = filter(grid, rank(dist) <= 10)
  )

ニュートン法による最良モデルの探索
モデルとデータセットとの距離を定義する関数とモデルのパラメータを変更して距離を最小にできるアルゴリズムがあれば、最良のモデルが求まる。

best <- optim(c(0,0), measure_distance, data = sim1)  
best$par
## [1] 4.222248 2.051204
ggplot(sim1, aes(x,y)) +
  geom_point(size =2, color = "grey30") +
  geom_abline(intercept = best$par[1], slope = best$par[2])

線形モデルの適合専用に設計されたlm()はoptim()を使わないで、線形モデルの数学的構造を活用する。lm()はフォーミュラ（formula）という特別な形式でモデルファミリーを指定する。y ~ xはy = a_1 + a_2 * xと同等。

sim1_mod <- lm(y ~ x, data = sim1)  ## 
coef(sim1_mod)
## (Intercept)           x 
##    4.220822    2.051533

5.2 モデルの可視化

5.2.1 予測

予測はモデルが補足したパターンを示唆。 data_grid()はデータのある領域を等間隔に覆う格子を生成。modeler::add_prediction()は引数にデータフレームとモデルを取り、モデルの予測をデータフレームの新たな列に追加する。

grid <- sim1 %>%
  data_grid(x)
grid
## # A tibble: 10 × 1
##        x
##    <int>
##  1     1
##  2     2
##  3     3
##  4     4
##  5     5
##  6     6
##  7     7
##  8     8
##  9     9
## 10    10
#予測を追加、列predが追加される
grid <- grid %>%
  add_predictions(sim1_mod)
grid
## # A tibble: 10 × 2
##        x  pred
##    <int> <dbl>
##  1     1  6.27
##  2     2  8.32
##  3     3 10.4 
##  4     4 12.4 
##  5     5 14.5 
##  6     6 16.5 
##  7     7 18.6 
##  8     8 20.6 
##  9     9 22.7 
## 10    10 24.7
ggplot(sim1, aes(x)) +   ##予測をプロット
  geom_point(aes(y=y)) +
  geom_line(
    aes(y = pred),
    data= grid,
    colour = "red",
    linewidth = 1
  )

5.3 残差

残差とは観測値と予測値の距離。残差はモデルが見逃しているパターンを伝える。modeler::add_residuals()は引数に元のデータセットとモデルを取り、データフレームに残差を追加。残差の平均は常に0になる。

sim1 <- sim1 %>%
  add_residuals(sim1_mod)  ##残差 residが追加される
sim1
## # A tibble: 30 × 3
##        x     y    resid
##    <int> <dbl>    <dbl>
##  1     1  4.20 -2.07   
##  2     1  7.51  1.24   
##  3     1  2.13 -4.15   
##  4     2  8.99  0.665  
##  5     2 10.2   1.92   
##  6     2 11.3   2.97   
##  7     3  7.36 -3.02   
##  8     3 10.5   0.130  
##  9     3 10.5   0.136  
## 10     4 12.4   0.00763
## # ℹ 20 more rows
# 度数分布多角形を描いて残差の分布の広がりを理解する
ggplot(sim1, aes(resid)) +
  geom_freqpoly(binwidth = 0.5)

元の予測の代わりに残差を使ってプロットすることも多い。残差がランダムなノイズに見えるなら、モデルはデータセットのパターンをしっかり補足したと言える。

ggplot(sim1, aes(x, resid)) +
  geom_ref_line(h = 0) + ## 参照となる線を引く
  geom_point()

5.4 フォーミュラとモデルファミリー

model_matrix()はデータフレームとフォーミュラを引数にし、モデル方程式を表すtibbleを返す。出力の各列にはモデルの係数が関連づけられ、関数は常にy = a_1 * out1 + a_2 * out_2。

df <- tribble(
  ~y, ~x1, ~x2,
  4, 2, 5,
  5, 1, 6
)
# Rが切片をモデルに追加する時は、デフォルトで1ばかりの列を追加する
model_matrix(df, y ~ x1)
## # A tibble: 2 × 2
##   `(Intercept)`    x1
##           <dbl> <dbl>
## 1             1     2
## 2             1     1
# 1ばかりの列を追加したくなければ、-1で明示的に取り消す
model_matrix(df, y ~ x1 - 1)
## # A tibble: 2 × 1
##      x1
##   <dbl>
## 1     2
## 2     1
model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 × 3
##   `(Intercept)`    x1    x2
##           <dbl> <dbl> <dbl>
## 1             1     2     5
## 2             1     1     6

5.5 カテゴリ変数

model_matrix()はカテゴリ変数を(0, 1)の数値に変換できる。下記の例ではカテゴリ変数sexをsex_maleに変換。sex_maleはsexがmaleならば1、その他ならば0。

df <- tribble(
  ~ sex, ~ response,
  "male", 1,
  "female", 2,
  "male", 1
)
model_matrix(df, response ~ sex)
## # A tibble: 3 × 2
##   `(Intercept)` sexmale
##           <dbl>   <dbl>
## 1             1       1
## 2             1       0
## 3             1       1

sim2はカテゴリ変数x(a, b, c, d)を含むサンプルデータ。カテゴリ変数xのモデルは、各カテゴリで平均値を予測する。平均が偏差の二乗平均平方根距離（RMSD）を最小化するから。

data(sim2)
ggplot(sim2) +
  geom_point(aes(x,y))

mod2 <- lm(y ~ x, data = sim2) ## モデルを適合させて予測を生成
grid <- sim2 %>%
  data_grid(x) %>% 
  add_predictions(mod2)
grid
## # A tibble: 4 × 2
##   x      pred
##   <chr> <dbl>
## 1 a      1.15
## 2 b      8.12
## 3 c      6.13
## 4 d      1.91
ggplot(sim2, aes(x)) +
  geom_point(aes(y = y)) +
  geom_point(
    data = grid,
    aes(y = pred),
    color = "red",
    size =4
  )

5.6 交互作用（連続とカテゴリ）

sim3にはカテゴリ予測子x2と連続予測子x1が含まれる。

data(sim3)
ggplot(sim3, aes(x1, y)) +
  geom_point(aes(color = x2))

mod1：変数を+で追加するとモデルは効果を他の全てとは独立に推定。mod2：*で追加すると独立作用に加えて交互作用を考慮する。すなわちy ~ x1 * x2はy = a_0 + a_1 * x1 + a2_ * x2 + a_12 * x1 * x2、に翻訳される。

mod1 <- lm(y ~ x1 + x2, data = sim3)  ## 変数を+で追加するとモデルは効果を他の全てとは独立に推定
mod2 <- lm(y ~ x1 * x2, data = sim3)  ## *で追加すると独立作用に加えて交互作用を考慮する

2つの予測子があるのでdata_grid()に両方の変数を指定する必要がある、x1とx2の全ての一意な値を探し出し、全ての組み合わせを生成。両方のモデルから同時に予測を生成するには、spread_prediction()で各予測に対応する新たな列を作る。

grid <- sim3 %>% 
  data_grid(x1, x2) %>% 
  spread_predictions(mod1, mod2)
grid
## # A tibble: 40 × 4
##       x1 x2     mod1  mod2
##    <int> <fct> <dbl> <dbl>
##  1     1 a      1.67  1.21
##  2     1 b      4.56  7.52
##  3     1 c      6.48  5.71
##  4     1 d      4.03  2.32
##  5     2 a      1.48  1.12
##  6     2 b      4.37  6.66
##  7     2 c      6.28  5.68
##  8     2 d      3.84  2.50
##  9     3 a      1.28  1.02
## 10     3 b      4.17  5.81
## # ℹ 30 more rows

各予測を1列predにまとめて出力するには、行に追加するgather_prediction()を使用。

grid <- sim3 %>% 
  data_grid(x1, x2) %>% 
  gather_predictions(mod1, mod2)
grid
## # A tibble: 80 × 4
##    model    x1 x2     pred
##    <chr> <int> <fct> <dbl>
##  1 mod1      1 a      1.67
##  2 mod1      1 b      4.56
##  3 mod1      1 c      6.48
##  4 mod1      1 d      4.03
##  5 mod1      2 a      1.48
##  6 mod1      2 b      4.37
##  7 mod1      2 c      6.28
##  8 mod1      2 d      3.84
##  9 mod1      3 a      1.28
## 10 mod1      3 b      4.17
## # ℹ 70 more rows

両方のモデルの結果をfacet_wrap()で1つのプロットにまとめる。2つの変数の独立作用のみを考慮したmod1ではどの線の傾きも等しい。交互作用を考慮したmod2では各線で傾きも切片も異なる。

ggplot(sim3, aes(x1, y, color = x2)) +
  geom_point() +
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model)

残差を調べて、どちらのモデルが適切かを検証。数学的手法ではなく、モデルが対象としているパターンを補足したかどうかの定性的評価に着目。mod2の残差では自明なパターンがほとんど無いのに対し、mod1の残差では、bである種のパターンが見逃されている。

sim3 <- sim3 %>% 
  gather_residuals(mod1, mod2)
ggplot(sim3, aes(x1, resid, color =x2)) +
  geom_point() +
  facet_grid(model~x2)

5.7 交互作用（２つの連続変数）

sim4は2つの連続変数x1, x2を含むサンプルデータ。交互作用を考慮しない場合とする場合でモデルを構築。

data(sim4)  ## 2つの連続変数x1, x2を含む
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)

grid <- sim4 %>% 
  data_grid(
    x1 = seq_range(x1, 5), ## x1の一意の値すべてを使う代わりに最小値と最大値の間の5つの格子値を等間隔に使用
    x2 = seq_range(x2, 5) 
  ) %>% 
  gather_predictions(mod1, mod2)
grid
## # A tibble: 50 × 4
##    model    x1    x2   pred
##    <chr> <dbl> <dbl>  <dbl>
##  1 mod1   -1    -1    0.996
##  2 mod1   -1    -0.5 -0.395
##  3 mod1   -1     0   -1.79 
##  4 mod1   -1     0.5 -3.18 
##  5 mod1   -1     1   -4.57 
##  6 mod1   -0.5  -1    1.91 
##  7 mod1   -0.5  -0.5  0.516
##  8 mod1   -0.5   0   -0.875
##  9 mod1   -0.5   0.5 -2.27 
## 10 mod1   -0.5   1   -3.66 
## # ℹ 40 more rows

モデルを可視化、2つの連続予測子があるのでモデルを3次元表面で考える。

ggplot(grid, aes(x1, x2)) +
  geom_tile(aes(fill = pred)) +
  facet_wrap( ~ model)

表面を上から見るのではなく、複数のスライスを表示して横から見る。

ggplot(grid, aes(x1, pred, colour = x2, group = x2)) + 
  geom_line() +
  facet_wrap(~ model)

ggplot(grid, aes(x2, pred, colour = x1, group = x1)) + 
  geom_line() +
  facet_wrap(~ model)

5.8 変換

log(y) ~ sqrt(x1) + sqrt(x2)はlog(y) = a_1 + a_2 * sqrt(x1) + a_3 * sqrt(x2)、に変換される。
変換が+,*,^,-を含むならI()でラップして、Rがモデル指定の一部と間違えないようにする。
例： y ~ x + I(x^2)はy = a_1 + a_2*x + a_3 * x^2、と翻訳される。
I()を忘れて、y ~ x^2 + x、と指定するとRはy ~ x*x + xを計算する、x*xはxと自分自身の交互作用を意味するのでxと同じになり, y = a_1 + a_2*xを指定することになる。 モデルが何をしているのかは、model_matrix()で方程式lm()が適合しているのかは正確に何であるかを確認できる。

df <- tribble(
  ~y, ~x,
  1,  1,
  2,  2, 
  3,  3
)
model_matrix(df, y ~ x^2 + x)
## # A tibble: 3 × 2
##   `(Intercept)`     x
##           <dbl> <dbl>
## 1             1     1
## 2             1     2
## 3             1     3
model_matrix(df, y ~ I(x^2) + x)
## # A tibble: 3 × 3
##   `(Intercept)` `I(x^2)`     x
##           <dbl>    <dbl> <dbl>
## 1             1        1     1
## 2             1        4     2
## 3             1        9     3

poly()で多項式に適合させることができる。ただしデータの範囲外では、多項式が急激に正または負の無限大になるため、自然スプラインspline::ns()を使う方が安全。

model_matrix(df, y ~ poly(x, degree =2))
## # A tibble: 3 × 3
##   `(Intercept)` `poly(x, degree = 2)1` `poly(x, degree = 2)2`
##           <dbl>                  <dbl>                  <dbl>
## 1             1              -7.07e- 1                  0.408
## 2             1              -9.07e-17                 -0.816
## 3             1               7.07e- 1                  0.408
library(splines)
model_matrix(df, y ~ ns(x, df = 2)) # dfはdegree of freedom
## # A tibble: 3 × 3
##   `(Intercept)` `ns(x, df = 2)1` `ns(x, df = 2)2`
##           <dbl>            <dbl>            <dbl>
## 1             1            0                0    
## 2             1            0.566           -0.211
## 3             1            0.344            0.771

非線形関数の近似。

sim5 <- tibble(
  x = seq(0, 3.5*pi, length = 50),
  y = 4 * sin(x) + rnorm(length(x))
)
ggplot(sim5, aes(x,y)) + geom_point()

# 以下の5つのモデルを適合させる
mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)

grid <- sim5 %>% 
  data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>%  ## expand = 0.1は範囲を10%拡張する
  gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")

ggplot(sim5, aes(x, y)) + 
  geom_point() +
  geom_line(data = grid, colour = "red") +
  facet_wrap(~ model)

5.9 欠損値

欠損値は変数間の関係に何の情報ももたらさないので、モデル関数は欠損値を含む行を削除する。

5.10 他のモデルファミリー

線形モデルでは残差が正規分布であることも仮定する。
一般化線形モデルglm()では目的変数の残差が正規分布に従わないデータに適用する、目的変数を変換して解析。

6 データ解析の実例

6.1 データの正規性の検証

現実のデータを統計解析する際は、まずデータが正規分布に従っているか否か、すなわち正規性があるかどうかを確認する。
まずデータをプロットして分布の状態を目視で確認するのと同時に、正規性の検定を行う。
データの正規性の検定にはシャピロ・ウィルク検定がよく用いられる。Rではデフォルトでインストールされている関数shapiro.test()で検定することができる。この検定の帰無仮説は「データxが正規分布している」であり、計算されたp値が0.05以上ならば帰無仮説を棄却できず、従ってデータxは正規分布していると結論される。

データxが正規分布する例：

x <-rnorm(50, mean = 5, sd = 3)
shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.9824, p-value = 0.6568

データxが正規分布しない例：

x <- runif(50, min = 0, max = 100)
shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.94977, p-value = 0.03337

shapiro.test()に渡すデータにはNaNが含まれていても計算できるが、NaN以外の要素の個数が3〜5000という制限がある。そのため5000を超える場合はsample()による無作為抽出で対応する。

x <- rnorm(6000, mean = 10, sd = 3)
sample_5000 <- sample(x, 5000, replace = FALSE) 
shapiro.test(sample_5000)
## 
##  Shapiro-Wilk normality test
## 
## data:  sample_5000
## W = 0.99962, p-value = 0.4695

shapiro.test()の入力はベクトルでありデータフレーム形式には対応していないため、データフレーム中のデータを渡すには自作関数を作る必要がある。実例はこちらのページ。

6.2 パラメトリック検定とノンパラメトリック検定

2群のデータを比較する場合、データに正規性がある場合はパラメトリック検定、正規性が無い場合はノンパラメトリック検定で統計検定する。

パラメトリック検定

  • Studentのt検定（2群比較で等分散を仮定）
  • Welchのt検定（2群比較で等分散を仮定しない）

ノンパラメトリック検定

  • Mann-Whitneyのt検定（2群比較で等分散を仮定しない）

6.3 多群データの比較

3群以上のデータを比較する場合、データの正規性の検証後にさらに一元配置分散分析により各群の平均値に差があるか検証する必要がある。有意な差が検出された場合は、事後検定で比較する。

正規性がある場合

  • ANOVA（パラメトリックな一元配置分散分析）
  • Tukey-Kramer検定（全ての群の組み合わせを比較）
  • Dunnett検定（対照群とそれ以外を比較）

正規性が無い場合

  • Kruskal-Wallis検定（ノンパラメトリックな一元配置分散分析）
  • Steel-Dwass検定（全ての群の組み合わせを比較）
  • Steel検定（対照群とそれ以外を比較）

各種検定の計算方法については、下記の実例集を参考にする。